Actual source code: plexgmsh.c
petsc-3.13.6 2020-09-29
1: #define PETSCDM_DLL
2: #include <petsc/private/dmpleximpl.h>
3: #include <petsc/private/hashmapi.h>
5: typedef struct {
6: PetscViewer viewer;
7: int fileFormat;
8: int dataSize;
9: PetscBool binary;
10: PetscBool byteSwap;
11: size_t wlen;
12: void *wbuf;
13: size_t slen;
14: void *sbuf;
15: } GmshFile;
17: static PetscErrorCode GmshBufferGet(GmshFile *gmsh, size_t count, size_t eltsize, void *buf)
18: {
19: size_t size = count * eltsize;
23: if (gmsh->wlen < size) {
24: PetscFree(gmsh->wbuf);
25: PetscMalloc(size, &gmsh->wbuf);
26: gmsh->wlen = size;
27: }
28: *(void**)buf = size ? gmsh->wbuf : NULL;
29: return(0);
30: }
32: static PetscErrorCode GmshBufferSizeGet(GmshFile *gmsh, size_t count, void *buf)
33: {
34: size_t dataSize = (size_t)gmsh->dataSize;
35: size_t size = count * dataSize;
39: if (gmsh->slen < size) {
40: PetscFree(gmsh->sbuf);
41: PetscMalloc(size, &gmsh->sbuf);
42: gmsh->slen = size;
43: }
44: *(void**)buf = size ? gmsh->sbuf : NULL;
45: return(0);
46: }
48: static PetscErrorCode GmshRead(GmshFile *gmsh, void *buf, PetscInt count, PetscDataType dtype)
49: {
52: PetscViewerRead(gmsh->viewer, buf, count, NULL, dtype);
53: if (gmsh->byteSwap) {PetscByteSwap(buf, dtype, count);}
54: return(0);
55: }
57: static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscInt count)
58: {
61: PetscViewerRead(gmsh->viewer, buf, count, NULL, PETSC_STRING);
62: return(0);
63: }
65: static PetscErrorCode GmshMatch(PETSC_UNUSED GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN], PetscBool *match)
66: {
69: PetscStrcmp(line, Section, match);
70: return(0);
71: }
73: static PetscErrorCode GmshExpect(GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN])
74: {
75: PetscBool match;
79: GmshMatch(gmsh, Section, line, &match);
80: if (!match) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file, expecting %s",Section);
81: return(0);
82: }
84: static PetscErrorCode GmshReadSection(GmshFile *gmsh, char line[PETSC_MAX_PATH_LEN])
85: {
86: PetscBool match;
90: while (PETSC_TRUE) {
91: GmshReadString(gmsh, line, 1);
92: GmshMatch(gmsh, "$Comments", line, &match);
93: if (!match) break;
94: while (PETSC_TRUE) {
95: GmshReadString(gmsh, line, 1);
96: GmshMatch(gmsh, "$EndComments", line, &match);
97: if (match) break;
98: }
99: }
100: return(0);
101: }
103: static PetscErrorCode GmshReadEndSection(GmshFile *gmsh, const char EndSection[], char line[PETSC_MAX_PATH_LEN])
104: {
107: GmshReadString(gmsh, line, 1);
108: GmshExpect(gmsh, EndSection, line);
109: return(0);
110: }
112: static PetscErrorCode GmshReadSize(GmshFile *gmsh, PetscInt *buf, PetscInt count)
113: {
114: PetscInt i;
115: size_t dataSize = (size_t)gmsh->dataSize;
119: if (dataSize == sizeof(PetscInt)) {
120: GmshRead(gmsh, buf, count, PETSC_INT);
121: } else if (dataSize == sizeof(int)) {
122: int *ibuf = NULL;
123: GmshBufferSizeGet(gmsh, count, &ibuf);
124: GmshRead(gmsh, ibuf, count, PETSC_ENUM);
125: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
126: } else if (dataSize == sizeof(long)) {
127: long *ibuf = NULL;
128: GmshBufferSizeGet(gmsh, count, &ibuf);
129: GmshRead(gmsh, ibuf, count, PETSC_LONG);
130: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
131: } else if (dataSize == sizeof(PetscInt64)) {
132: PetscInt64 *ibuf = NULL;
133: GmshBufferSizeGet(gmsh, count, &ibuf);
134: GmshRead(gmsh, ibuf, count, PETSC_INT64);
135: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
136: }
137: return(0);
138: }
140: static PetscErrorCode GmshReadInt(GmshFile *gmsh, int *buf, PetscInt count)
141: {
144: GmshRead(gmsh, buf, count, PETSC_ENUM);
145: return(0);
146: }
148: static PetscErrorCode GmshReadDouble(GmshFile *gmsh, double *buf, PetscInt count)
149: {
152: GmshRead(gmsh, buf, count, PETSC_DOUBLE);
153: return(0);
154: }
156: typedef struct {
157: PetscInt id; /* Entity number */
158: PetscInt dim; /* Entity dimension */
159: double bbox[6]; /* Bounding box */
160: PetscInt numTags; /* Size of tag array */
161: int tags[4]; /* Tag array */
162: } GmshEntity;
164: typedef struct {
165: GmshEntity *entity[4];
166: PetscHMapI entityMap[4];
167: } GmshEntities;
169: static PetscErrorCode GmshEntitiesCreate(PetscInt count[4], GmshEntities **entities)
170: {
171: PetscInt dim;
175: PetscNew(entities);
176: for (dim = 0; dim < 4; ++dim) {
177: PetscCalloc1(count[dim], &(*entities)->entity[dim]);
178: PetscHMapICreate(&(*entities)->entityMap[dim]);
179: }
180: return(0);
181: }
183: static PetscErrorCode GmshEntitiesAdd(GmshEntities *entities, PetscInt index, PetscInt dim, PetscInt eid, GmshEntity** entity)
184: {
187: PetscHMapISet(entities->entityMap[dim], eid, index);
188: entities->entity[dim][index].dim = dim;
189: entities->entity[dim][index].id = eid;
190: if (entity) *entity = &entities->entity[dim][index];
191: return(0);
192: }
194: static PetscErrorCode GmshEntitiesGet(GmshEntities *entities, PetscInt dim, PetscInt eid, GmshEntity** entity)
195: {
196: PetscInt index;
200: PetscHMapIGet(entities->entityMap[dim], eid, &index);
201: *entity = &entities->entity[dim][index];
202: return(0);
203: }
205: static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities)
206: {
207: PetscInt dim;
211: if (!*entities) return(0);
212: for (dim = 0; dim < 4; ++dim) {
213: PetscFree((*entities)->entity[dim]);
214: PetscHMapIDestroy(&(*entities)->entityMap[dim]);
215: }
216: PetscFree((*entities));
217: return(0);
218: }
220: typedef struct {
221: PetscInt id; /* Entity number */
222: PetscInt dim; /* Entity dimension */
223: PetscInt cellType; /* Cell type */
224: PetscInt numNodes; /* Size of node array */
225: PetscInt nodes[8]; /* Node array */
226: PetscInt numTags; /* Size of tag array */
227: int tags[4]; /* Tag array */
228: } GmshElement;
230: static PetscErrorCode DMPlexCreateGmsh_ReadNodes_v22(GmshFile *gmsh, int shift, PetscInt *numVertices, double **gmsh_nodes)
231: {
232: PetscViewer viewer = gmsh->viewer;
233: PetscBool byteSwap = gmsh->byteSwap;
234: char line[PETSC_MAX_PATH_LEN];
235: int v, num, nid, snum;
236: double *coordinates;
240: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
241: snum = sscanf(line, "%d", &num);
242: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
243: PetscMalloc1(num*3, &coordinates);
244: *numVertices = num;
245: *gmsh_nodes = coordinates;
246: for (v = 0; v < num; ++v) {
247: double *xyz = coordinates + v*3;
248: PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
249: if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
250: if (nid != v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %d", nid, v+shift);
251: PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
252: if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
253: }
254: return(0);
255: }
257: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
258: file contents multiple times to figure out the true number of cells and facets
259: in the given mesh. To make this more efficient we read the file contents only
260: once and store them in memory, while determining the true number of cells. */
261: static PetscErrorCode DMPlexCreateGmsh_ReadElements_v22(GmshFile* gmsh, PETSC_UNUSED int shift, PetscInt *numCells, GmshElement **gmsh_elems)
262: {
263: PetscViewer viewer = gmsh->viewer;
264: PetscBool binary = gmsh->binary;
265: PetscBool byteSwap = gmsh->byteSwap;
266: char line[PETSC_MAX_PATH_LEN];
267: GmshElement *elements;
268: int i, c, p, num, ibuf[1+4+512], snum;
269: int cellType, dim, numNodes, numNodesIgnore, numElem, numTags;
273: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
274: snum = sscanf(line, "%d", &num);
275: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
276: PetscMalloc1(num, &elements);
277: *numCells = num;
278: *gmsh_elems = elements;
279: for (c = 0; c < num;) {
280: PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
281: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 3);}
282: if (binary) {
283: cellType = ibuf[0];
284: numElem = ibuf[1];
285: numTags = ibuf[2];
286: } else {
287: elements[c].id = ibuf[0];
288: cellType = ibuf[1];
289: numTags = ibuf[2];
290: numElem = 1;
291: }
292: /* http://gmsh.info/doc/texinfo/gmsh.html#MSH-ASCII-file-format */
293: numNodesIgnore = 0;
294: switch (cellType) {
295: case 1: /* 2-node line */
296: dim = 1;
297: numNodes = 2;
298: break;
299: case 2: /* 3-node triangle */
300: dim = 2;
301: numNodes = 3;
302: break;
303: case 3: /* 4-node quadrangle */
304: dim = 2;
305: numNodes = 4;
306: break;
307: case 4: /* 4-node tetrahedron */
308: dim = 3;
309: numNodes = 4;
310: break;
311: case 5: /* 8-node hexahedron */
312: dim = 3;
313: numNodes = 8;
314: break;
315: case 6: /* 6-node wedge */
316: dim = 3;
317: numNodes = 6;
318: break;
319: case 8: /* 3-node 2nd order line */
320: dim = 1;
321: numNodes = 2;
322: numNodesIgnore = 1;
323: break;
324: case 9: /* 6-node 2nd order triangle */
325: dim = 2;
326: numNodes = 3;
327: numNodesIgnore = 3;
328: break;
329: case 10: /* 9-node 2nd order quadrangle */
330: dim = 2;
331: numNodes = 4;
332: numNodesIgnore = 5;
333: break;
334: case 11: /* 10-node 2nd order tetrahedron */
335: dim = 3;
336: numNodes = 4;
337: numNodesIgnore = 6;
338: break;
339: case 12: /* 27-node 2nd order hexhedron */
340: dim = 3;
341: numNodes = 8;
342: numNodesIgnore = 19;
343: break;
344: case 13: /* 18-node 2nd wedge */
345: dim = 3;
346: numNodes = 6;
347: numNodesIgnore = 12;
348: break;
349: case 15: /* 1-node vertex */
350: dim = 0;
351: numNodes = 1;
352: break;
353: case 7: /* 5-node pyramid */
354: case 14: /* 14-node 2nd order pyramid */
355: default:
356: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
357: }
358: if (binary) {
359: const int nint = 1 + numTags + numNodes + numNodesIgnore;
360: /* Loop over element blocks */
361: for (i = 0; i < numElem; ++i, ++c) {
362: PetscViewerRead(viewer, ibuf, nint, NULL, PETSC_ENUM);
363: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, nint);}
364: elements[c].dim = dim;
365: elements[c].numNodes = numNodes;
366: elements[c].numTags = numTags;
367: elements[c].id = ibuf[0];
368: elements[c].cellType = cellType;
369: for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p];
370: for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
371: }
372: } else {
373: const int nint = numTags + numNodes + numNodesIgnore;
374: elements[c].dim = dim;
375: elements[c].numNodes = numNodes;
376: elements[c].numTags = numTags;
377: elements[c].cellType = cellType;
378: PetscViewerRead(viewer, ibuf, nint, NULL, PETSC_ENUM);
379: for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[p];
380: for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[numTags + p];
381: c++;
382: }
383: }
384: return(0);
385: }
387: /*
388: $Entities
389: numPoints(unsigned long) numCurves(unsigned long)
390: numSurfaces(unsigned long) numVolumes(unsigned long)
391: // points
392: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
393: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
394: numPhysicals(unsigned long) phyisicalTag[...](int)
395: ...
396: // curves
397: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
398: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
399: numPhysicals(unsigned long) physicalTag[...](int)
400: numBREPVert(unsigned long) tagBREPVert[...](int)
401: ...
402: // surfaces
403: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
404: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
405: numPhysicals(unsigned long) physicalTag[...](int)
406: numBREPCurve(unsigned long) tagBREPCurve[...](int)
407: ...
408: // volumes
409: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
410: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
411: numPhysicals(unsigned long) physicalTag[...](int)
412: numBREPSurfaces(unsigned long) tagBREPSurfaces[...](int)
413: ...
414: $EndEntities
415: */
416: static PetscErrorCode DMPlexCreateGmsh_ReadEntities_v40(GmshFile *gmsh, GmshEntities **entities)
417: {
418: PetscViewer viewer = gmsh->viewer;
419: PetscBool byteSwap = gmsh->byteSwap;
420: long index, num, lbuf[4];
421: int dim, eid, numTags, *ibuf, t;
422: PetscInt count[4], i;
423: GmshEntity *entity = NULL;
427: PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG);
428: if (byteSwap) {PetscByteSwap(lbuf, PETSC_LONG, 4);}
429: for (i = 0; i < 4; ++i) count[i] = lbuf[i];
430: GmshEntitiesCreate(count, entities);
431: for (dim = 0; dim < 4; ++dim) {
432: for (index = 0; index < count[dim]; ++index) {
433: PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM);
434: if (byteSwap) {PetscByteSwap(&eid, PETSC_ENUM, 1);}
435: GmshEntitiesAdd(*entities, (PetscInt)index, dim, eid, &entity);
436: PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE);
437: if (byteSwap) {PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6);}
438: PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
439: if (byteSwap) {PetscByteSwap(&num, PETSC_LONG, 1);}
440: GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
441: PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
442: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, num);}
443: entity->numTags = numTags = (int) PetscMin(num, 4);
444: for (t = 0; t < numTags; ++t) entity->tags[t] = ibuf[t];
445: if (dim == 0) continue;
446: PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
447: if (byteSwap) {PetscByteSwap(&num, PETSC_LONG, 1);}
448: GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
449: PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
450: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, num);}
451: }
452: }
453: return(0);
454: }
456: /*
457: $Nodes
458: numEntityBlocks(unsigned long) numNodes(unsigned long)
459: tagEntity(int) dimEntity(int) typeNode(int) numNodes(unsigned long)
460: tag(int) x(double) y(double) z(double)
461: ...
462: ...
463: $EndNodes
464: */
465: static PetscErrorCode DMPlexCreateGmsh_ReadNodes_v40(GmshFile *gmsh, int shift, PetscInt *numVertices, double **gmsh_nodes)
466: {
467: PetscViewer viewer = gmsh->viewer;
468: PetscBool byteSwap = gmsh->byteSwap;
469: long block, node, v, numEntityBlocks, numTotalNodes, numNodes;
470: int info[3], nid;
471: double *coordinates;
475: PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
476: if (byteSwap) {PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);}
477: PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG);
478: if (byteSwap) {PetscByteSwap(&numTotalNodes, PETSC_LONG, 1);}
479: PetscMalloc1(numTotalNodes*3, &coordinates);
480: *numVertices = numTotalNodes;
481: *gmsh_nodes = coordinates;
482: for (v = 0, block = 0; block < numEntityBlocks; ++block) {
483: PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
484: PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG);
485: if (byteSwap) {PetscByteSwap(&numNodes, PETSC_LONG, 1);}
486: if (gmsh->binary) {
487: size_t nbytes = sizeof(int) + 3*sizeof(double);
488: char *cbuf = NULL; /* dummy value to prevent warning from compiler about possible unitilized value */
489: GmshBufferGet(gmsh, numNodes, nbytes, &cbuf);
490: PetscViewerRead(viewer, cbuf, numNodes*nbytes, NULL, PETSC_CHAR);
491: for (node = 0; node < numNodes; ++node, ++v) {
492: char *cnid = cbuf + node*nbytes, *cxyz = cnid + sizeof(int);
493: double *xyz = coordinates + v*3;
494: if (!PetscBinaryBigEndian()) {PetscByteSwap(cnid, PETSC_ENUM, 1);}
495: if (!PetscBinaryBigEndian()) {PetscByteSwap(cxyz, PETSC_DOUBLE, 3);}
496: PetscMemcpy(&nid, cnid, sizeof(int));
497: PetscMemcpy(xyz, cxyz, 3*sizeof(double));
498: if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
499: if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
500: if ((long)nid != v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %ld", nid, v+shift);
501: }
502: } else {
503: for (node = 0; node < numNodes; ++node, ++v) {
504: double *xyz = coordinates + v*3;
505: PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
506: if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
507: if ((long)nid != v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %ld", nid, v+shift);
508: PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
509: if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
510: }
511: }
512: }
513: return(0);
514: }
516: /*
517: $Elements
518: numEntityBlocks(unsigned long) numElements(unsigned long)
519: tagEntity(int) dimEntity(int) typeEle(int) numElements(unsigned long)
520: tag(int) numVert[...](int)
521: ...
522: ...
523: $EndElements
524: */
525: static PetscErrorCode DMPlexCreateGmsh_ReadElements_v40(GmshFile *gmsh, PETSC_UNUSED int shift, GmshEntities *entities, PetscInt *numCells, GmshElement **gmsh_elems)
526: {
527: PetscViewer viewer = gmsh->viewer;
528: PetscBool byteSwap = gmsh->byteSwap;
529: long c, block, numEntityBlocks, numTotalElements, elem, numElements;
530: int p, info[3], *ibuf = NULL;
531: int eid, dim, numTags, *tags, cellType, numNodes;
532: GmshEntity *entity = NULL;
533: GmshElement *elements;
537: PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
538: if (byteSwap) {PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);}
539: PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG);
540: if (byteSwap) {PetscByteSwap(&numTotalElements, PETSC_LONG, 1);}
541: PetscCalloc1(numTotalElements, &elements);
542: *numCells = numTotalElements;
543: *gmsh_elems = elements;
544: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
545: PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
546: if (byteSwap) {PetscByteSwap(info, PETSC_ENUM, 3);}
547: eid = info[0]; dim = info[1]; cellType = info[2];
548: GmshEntitiesGet(entities, dim, eid, &entity);
549: numTags = entity->numTags;
550: tags = entity->tags;
551: switch (cellType) {
552: case 1: /* 2-node line */
553: numNodes = 2;
554: break;
555: case 2: /* 3-node triangle */
556: numNodes = 3;
557: break;
558: case 3: /* 4-node quadrangle */
559: numNodes = 4;
560: break;
561: case 4: /* 4-node tetrahedron */
562: numNodes = 4;
563: break;
564: case 5: /* 8-node hexahedron */
565: numNodes = 8;
566: break;
567: case 6: /* 6-node wedge */
568: numNodes = 6;
569: break;
570: case 15: /* 1-node vertex */
571: numNodes = 1;
572: break;
573: default:
574: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
575: }
576: PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG);
577: if (byteSwap) {PetscByteSwap(&numElements, PETSC_LONG, 1);}
578: GmshBufferGet(gmsh, (1+numNodes)*numElements, sizeof(int), &ibuf);
579: PetscViewerRead(viewer, ibuf, (1+numNodes)*numElements, NULL, PETSC_ENUM);
580: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, (1+numNodes)*numElements);}
581: for (elem = 0; elem < numElements; ++elem, ++c) {
582: int *id = ibuf + elem*(1+numNodes), *nodes = id + 1;
583: GmshElement *element = elements + c;
584: element->dim = dim;
585: element->cellType = cellType;
586: element->numNodes = numNodes;
587: element->numTags = numTags;
588: element->id = id[0];
589: for (p = 0; p < numNodes; p++) element->nodes[p] = nodes[p];
590: for (p = 0; p < numTags; p++) element->tags[p] = tags[p];
591: }
592: }
593: return(0);
594: }
596: static PetscErrorCode DMPlexCreateGmsh_ReadPeriodic_v40(GmshFile *gmsh, int shift, PetscInt slaveMap[], PetscBT bt)
597: {
598: PetscViewer viewer = gmsh->viewer;
599: int fileFormat = gmsh->fileFormat;
600: PetscBool binary = gmsh->binary;
601: PetscBool byteSwap = gmsh->byteSwap;
602: int numPeriodic, snum, i;
603: char line[PETSC_MAX_PATH_LEN];
607: if (fileFormat == 22 || !binary) {
608: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
609: snum = sscanf(line, "%d", &numPeriodic);
610: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
611: } else {
612: PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM);
613: if (byteSwap) {PetscByteSwap(&numPeriodic, PETSC_ENUM, 1);}
614: }
615: for (i = 0; i < numPeriodic; i++) {
616: int ibuf[3], slaveDim = -1, slaveTag = -1, masterTag = -1, slaveNode, masterNode;
617: long j, nNodes;
618: double affine[16];
620: if (fileFormat == 22 || !binary) {
621: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
622: snum = sscanf(line, "%d %d %d", &slaveDim, &slaveTag, &masterTag);
623: if (snum != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
624: } else {
625: PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
626: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 3);}
627: slaveDim = ibuf[0]; slaveTag = ibuf[1]; masterTag = ibuf[2];
628: }
629: (void)slaveDim; (void)slaveTag; (void)masterTag; /* unused */
631: if (fileFormat == 22 || !binary) {
632: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
633: snum = sscanf(line, "%ld", &nNodes);
634: if (snum != 1) { /* discard transformation and try again */
635: PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
636: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
637: snum = sscanf(line, "%ld", &nNodes);
638: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
639: }
640: } else {
641: PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
642: if (byteSwap) {PetscByteSwap(&nNodes, PETSC_LONG, 1);}
643: if (nNodes == -1) { /* discard transformation and try again */
644: PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE);
645: PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
646: if (byteSwap) {PetscByteSwap(&nNodes, PETSC_LONG, 1);}
647: }
648: }
650: for (j = 0; j < nNodes; j++) {
651: if (fileFormat == 22 || !binary) {
652: PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
653: snum = sscanf(line, "%d %d", &slaveNode, &masterNode);
654: if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
655: } else {
656: PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM);
657: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 2);}
658: slaveNode = ibuf[0]; masterNode = ibuf[1];
659: }
660: slaveMap[slaveNode - shift] = masterNode - shift;
661: PetscBTSet(bt, slaveNode - shift);
662: PetscBTSet(bt, masterNode - shift);
663: }
664: }
665: return(0);
666: }
668: /*
669: $Entities
670: numPoints(size_t) numCurves(size_t)
671: numSurfaces(size_t) numVolumes(size_t)
672: pointTag(int) X(double) Y(double) Z(double)
673: numPhysicalTags(size_t) physicalTag(int) ...
674: ...
675: curveTag(int) minX(double) minY(double) minZ(double)
676: maxX(double) maxY(double) maxZ(double)
677: numPhysicalTags(size_t) physicalTag(int) ...
678: numBoundingPoints(size_t) pointTag(int) ...
679: ...
680: surfaceTag(int) minX(double) minY(double) minZ(double)
681: maxX(double) maxY(double) maxZ(double)
682: numPhysicalTags(size_t) physicalTag(int) ...
683: numBoundingCurves(size_t) curveTag(int) ...
684: ...
685: volumeTag(int) minX(double) minY(double) minZ(double)
686: maxX(double) maxY(double) maxZ(double)
687: numPhysicalTags(size_t) physicalTag(int) ...
688: numBoundngSurfaces(size_t) surfaceTag(int) ...
689: ...
690: $EndEntities
691: */
692: static PetscErrorCode DMPlexCreateGmsh_ReadEntities_v41(GmshFile *gmsh, GmshEntities **entities)
693: {
694: PetscInt count[4], index, numTags, i;
695: int dim, eid, *tags = NULL;
696: GmshEntity *entity = NULL;
700: GmshReadSize(gmsh, count, 4);
701: GmshEntitiesCreate(count, entities);
702: for (dim = 0; dim < 4; ++dim) {
703: for (index = 0; index < count[dim]; ++index) {
704: GmshReadInt(gmsh, &eid, 1);
705: GmshEntitiesAdd(*entities, (PetscInt)index, dim, eid, &entity);
706: GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6);
707: GmshReadSize(gmsh, &numTags, 1);
708: GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
709: GmshReadInt(gmsh, tags, numTags);
710: entity->numTags = PetscMin(numTags, 4);
711: for (i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i];
712: if (dim == 0) continue;
713: GmshReadSize(gmsh, &numTags, 1);
714: GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
715: GmshReadInt(gmsh, tags, numTags);
716: }
717: }
718: return(0);
719: }
721: /*
722: $Nodes
723: numEntityBlocks(size_t) numNodes(size_t)
724: minNodeTag(size_t) maxNodeTag(size_t)
725: entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesBlock(size_t)
726: nodeTag(size_t)
727: ...
728: x(double) y(double) z(double)
729: < u(double; if parametric and entityDim = 1 or entityDim = 2) >
730: < v(double; if parametric and entityDim = 2) >
731: ...
732: ...
733: $EndNodes
734: */
735: static PetscErrorCode DMPlexCreateGmsh_ReadNodes_v41(GmshFile *gmsh, int shift, PetscInt *numVertices, double **gmsh_nodes)
736: {
737: int info[3];
738: PetscInt sizes[4], numEntityBlocks, numNodes, numNodesBlock = 0, *nodeTag = NULL, block, node, i;
739: double *coordinates;
743: GmshReadSize(gmsh, sizes, 4);
744: numEntityBlocks = sizes[0]; numNodes = sizes[1];
745: PetscMalloc1(numNodes*3, &coordinates);
746: *numVertices = numNodes;
747: *gmsh_nodes = coordinates;
748: for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) {
749: GmshReadInt(gmsh, info, 3);
750: GmshReadSize(gmsh, &numNodesBlock, 1);
751: GmshBufferGet(gmsh, numNodesBlock, sizeof(PetscInt), &nodeTag);
752: GmshReadSize(gmsh, nodeTag, numNodesBlock);
753: for (i = 0; i < numNodesBlock; ++i) if (nodeTag[i] != node+i+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %d", nodeTag[i], node+i+shift);
754: if (info[2] != 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parametric coordinates not supported");
755: GmshReadDouble(gmsh, coordinates+node*3, numNodesBlock*3);
756: }
757: return(0);
758: }
760: /*
761: $Elements
762: numEntityBlocks(size_t) numElements(size_t)
763: minElementTag(size_t) maxElementTag(size_t)
764: entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t)
765: elementTag(size_t) nodeTag(size_t) ...
766: ...
767: ...
768: $EndElements
769: */
770: static PetscErrorCode DMPlexCreateGmsh_ReadElements_v41(GmshFile *gmsh, PETSC_UNUSED int shift, GmshEntities *entities, PetscInt *numCells, GmshElement **gmsh_elems)
771: {
772: int info[3], eid, dim, cellType, *tags;
773: PetscInt sizes[4], *ibuf = NULL, numEntityBlocks, numElements, numBlockElements, numNodes, numTags, block, elem, c, p;
774: GmshEntity *entity = NULL;
775: GmshElement *elements;
779: GmshReadSize(gmsh, sizes, 4);
780: numEntityBlocks = sizes[0]; numElements = sizes[1];
781: PetscCalloc1(numElements, &elements);
782: *numCells = numElements;
783: *gmsh_elems = elements;
784: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
785: GmshReadInt(gmsh, info, 3);
786: dim = info[0]; eid = info[1]; cellType = info[2];
787: GmshEntitiesGet(entities, dim, eid, &entity);
788: numTags = entity->numTags;
789: tags = entity->tags;
790: switch (cellType) {
791: case 1: /* 2-node line */
792: numNodes = 2;
793: break;
794: case 2: /* 3-node triangle */
795: numNodes = 3;
796: break;
797: case 3: /* 4-node quadrangle */
798: numNodes = 4;
799: break;
800: case 4: /* 4-node tetrahedron */
801: numNodes = 4;
802: break;
803: case 5: /* 8-node hexahedron */
804: numNodes = 8;
805: break;
806: case 6: /* 6-node wedge */
807: numNodes = 6;
808: break;
809: case 15: /* 1-node vertex */
810: numNodes = 1;
811: break;
812: default:
813: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
814: }
815: GmshReadSize(gmsh, &numBlockElements, 1);
816: GmshBufferGet(gmsh, (1+numNodes)*numBlockElements, sizeof(PetscInt), &ibuf);
817: GmshReadSize(gmsh, ibuf, (1+numNodes)*numBlockElements);
818: for (elem = 0; elem < numBlockElements; ++elem, ++c) {
819: GmshElement *element = elements + c;
820: PetscInt *id = ibuf + elem*(1+numNodes), *nodes = id + 1;
821: element->id = id[0];
822: element->dim = dim;
823: element->cellType = cellType;
824: element->numNodes = numNodes;
825: element->numTags = numTags;
826: for (p = 0; p < numNodes; p++) element->nodes[p] = nodes[p];
827: for (p = 0; p < numTags; p++) element->tags[p] = tags[p];
828: }
829: }
830: return(0);
831: }
833: /*
834: $Periodic
835: numPeriodicLinks(size_t)
836: entityDim(int) entityTag(int) entityTagMaster(int)
837: numAffine(size_t) value(double) ...
838: numCorrespondingNodes(size_t)
839: nodeTag(size_t) nodeTagMaster(size_t)
840: ...
841: ...
842: $EndPeriodic
843: */
844: static PetscErrorCode DMPlexCreateGmsh_ReadPeriodic_v41(GmshFile *gmsh, int shift, PetscInt slaveMap[], PetscBT bt)
845: {
846: int info[3];
847: PetscInt numPeriodicLinks, numAffine, numCorrespondingNodes, *nodeTags = NULL, link, node;
848: double dbuf[16];
852: GmshReadSize(gmsh, &numPeriodicLinks, 1);
853: for (link = 0; link < numPeriodicLinks; ++link) {
854: GmshReadInt(gmsh, info, 3);
855: GmshReadSize(gmsh, &numAffine, 1);
856: GmshReadDouble(gmsh, dbuf, numAffine);
857: GmshReadSize(gmsh, &numCorrespondingNodes, 1);
858: GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags);
859: GmshReadSize(gmsh, nodeTags, numCorrespondingNodes*2);
860: for (node = 0; node < numCorrespondingNodes; ++node) {
861: PetscInt slaveNode = nodeTags[node*2+0] - shift;
862: PetscInt masterNode = nodeTags[node*2+1] - shift;
863: slaveMap[slaveNode] = masterNode;
864: PetscBTSet(bt, slaveNode);
865: PetscBTSet(bt, masterNode);
866: }
867: }
868: return(0);
869: }
871: /*
872: $MeshFormat // same as MSH version 2
873: version(ASCII double; currently 4.1)
874: file-type(ASCII int; 0 for ASCII mode, 1 for binary mode)
875: data-size(ASCII int; sizeof(size_t))
876: < int with value one; only in binary mode, to detect endianness >
877: $EndMeshFormat
878: */
879: static PetscErrorCode DMPlexCreateGmsh_ReadMeshFormat(GmshFile *gmsh)
880: {
881: char line[PETSC_MAX_PATH_LEN];
882: int snum, fileType, fileFormat, dataSize, checkEndian;
883: float version;
887: GmshReadString(gmsh, line, 3);
888: snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
889: if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
890: if (version < 2.2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at least 2.2", (double)version);
891: if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f not supported", (double)version);
892: if (version > 4.1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at most 4.1", (double)version);
893: if (gmsh->binary && !fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Viewer is binary but Gmsh file is ASCII");
894: if (!gmsh->binary && fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Viewer is ASCII but Gmsh file is binary");
895: fileFormat = (int)(version*10); /* XXX Should use (int)roundf(version*10) ? */
896: if (fileFormat <= 40 && dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
897: if (fileFormat >= 41 && dataSize != sizeof(int) && dataSize != sizeof(PetscInt64)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
898: gmsh->fileFormat = fileFormat;
899: gmsh->dataSize = dataSize;
900: gmsh->byteSwap = PETSC_FALSE;
901: if (gmsh->binary) {
902: GmshReadInt(gmsh, &checkEndian, 1);
903: if (checkEndian != 1) {
904: PetscByteSwap(&checkEndian, PETSC_ENUM, 1);
905: if (checkEndian != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to detect endianness in Gmsh file header: %s", line);
906: gmsh->byteSwap = PETSC_TRUE;
907: }
908: }
909: return(0);
910: }
912: /*
913: PhysicalNames
914: numPhysicalNames(ASCII int)
915: dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max)
916: ...
917: $EndPhysicalNames
918: */
919: static PetscErrorCode DMPlexCreateGmsh_ReadPhysicalNames(GmshFile *gmsh)
920: {
921: char line[PETSC_MAX_PATH_LEN], name[128+2], *p, *q;
922: int snum, numRegions, region, dim, tag;
926: GmshReadString(gmsh, line, 1);
927: snum = sscanf(line, "%d", &numRegions);
928: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
929: for (region = 0; region < numRegions; ++region) {
930: GmshReadString(gmsh, line, 2);
931: snum = sscanf(line, "%d %d", &dim, &tag);
932: if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
933: GmshReadString(gmsh, line, -(PetscInt)sizeof(line));
934: PetscStrchr(line, '"', &p);
935: if (!p) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
936: PetscStrrchr(line, '"', &q);
937: if (q == p) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
938: PetscStrncpy(name, p+1, (size_t)(q-p-1));
939: }
940: return(0);
941: }
944: PETSC_STATIC_INLINE DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt ctGmsh)
945: {
946: switch (ctGmsh) {
947: case 15:
948: return DM_POLYTOPE_POINT;
949: case 1:
950: case 8:
951: return DM_POLYTOPE_SEGMENT;
952: case 2:
953: case 9:
954: return DM_POLYTOPE_TRIANGLE;
955: case 3:
956: case 10:
957: return DM_POLYTOPE_QUADRILATERAL;
958: case 4:
959: case 11:
960: return DM_POLYTOPE_TETRAHEDRON;
961: case 5:
962: case 12:
963: return DM_POLYTOPE_HEXAHEDRON;
964: case 6:
965: case 13:
966: return DM_POLYTOPE_TRI_PRISM;
967: case 7:
968: case 14:
969: return DM_POLYTOPE_UNKNOWN; /* Pyramid */
970: default:
971: return DM_POLYTOPE_UNKNOWN;
972: }
973: }
975: /*@C
976: DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file
978: + comm - The MPI communicator
979: . filename - Name of the Gmsh file
980: - interpolate - Create faces and edges in the mesh
982: Output Parameter:
983: . dm - The DM object representing the mesh
985: Level: beginner
987: .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
988: @*/
989: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
990: {
991: PetscViewer viewer;
992: PetscMPIInt rank;
993: int fileType;
994: PetscViewerType vtype;
995: PetscErrorCode ierr;
998: MPI_Comm_rank(comm, &rank);
1000: /* Determine Gmsh file type (ASCII or binary) from file header */
1001: if (!rank) {
1002: GmshFile gmsh_, *gmsh = &gmsh_;
1003: char line[PETSC_MAX_PATH_LEN];
1004: int snum;
1005: float version;
1007: PetscArrayzero(gmsh,1);
1008: PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer);
1009: PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII);
1010: PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ);
1011: PetscViewerFileSetName(gmsh->viewer, filename);
1012: /* Read only the first two lines of the Gmsh file */
1013: GmshReadSection(gmsh, line);
1014: GmshExpect(gmsh, "$MeshFormat", line);
1015: GmshReadString(gmsh, line, 2);
1016: snum = sscanf(line, "%f %d", &version, &fileType);
1017: if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
1018: if (version < 2.2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at least 2.2", (double)version);
1019: if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f not supported", (double)version);
1020: if (version > 4.1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at most 4.1", (double)version);
1021: PetscViewerDestroy(&gmsh->viewer);
1022: }
1023: MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);
1024: vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;
1026: /* Create appropriate viewer and build plex */
1027: PetscViewerCreate(comm, &viewer);
1028: PetscViewerSetType(viewer, vtype);
1029: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
1030: PetscViewerFileSetName(viewer, filename);
1031: DMPlexCreateGmsh(comm, viewer, interpolate, dm);
1032: PetscViewerDestroy(&viewer);
1033: return(0);
1034: }
1036: /*@
1037: DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer
1039: Collective
1041: Input Parameters:
1042: + comm - The MPI communicator
1043: . viewer - The Viewer associated with a Gmsh file
1044: - interpolate - Create faces and edges in the mesh
1046: Output Parameter:
1047: . dm - The DM object representing the mesh
1049: Note: http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
1051: Level: beginner
1053: .seealso: DMPLEX, DMCreate()
1054: @*/
1055: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1056: {
1057: PetscViewer parentviewer = NULL;
1058: double *coordsIn = NULL;
1059: GmshEntities *entities = NULL;
1060: GmshElement *gmsh_elem = NULL;
1061: PetscSection coordSection;
1062: Vec coordinates;
1063: PetscBT periodicV = NULL, periodicC = NULL;
1064: PetscScalar *coords;
1065: PetscInt dim = 0, embedDim = -1, coordSize, c, v, d, cell, *periodicMap = NULL, *periodicMapI = NULL, *hybridMap = NULL;
1066: PetscInt numVertices = 0, numCells = 0, trueNumCells = 0;
1067: int i, shift = 1;
1068: PetscMPIInt rank;
1069: PetscBool binary, zerobase = PETSC_FALSE, usemarker = PETSC_FALSE;
1070: PetscBool enable_hybrid = interpolate, periodic = PETSC_TRUE;
1071: PetscBool hasTetra = PETSC_FALSE;
1075: MPI_Comm_rank(comm, &rank);
1076: PetscObjectOptionsBegin((PetscObject)viewer);
1077: PetscOptionsHead(PetscOptionsObject,"DMPlex Gmsh options");
1078: PetscOptionsBool("-dm_plex_gmsh_hybrid", "Generate hybrid cell bounds", "DMPlexCreateGmsh", enable_hybrid, &enable_hybrid, NULL);
1079: PetscOptionsBool("-dm_plex_gmsh_periodic","Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL);
1080: PetscOptionsBool("-dm_plex_gmsh_use_marker", "Generate marker label", "DMPlexCreateGmsh", usemarker, &usemarker, NULL);
1081: PetscOptionsBool("-dm_plex_gmsh_zero_base", "Read Gmsh file with zero base indices", "DMPlexCreateGmsh", zerobase, &zerobase, NULL);
1082: PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", embedDim, &embedDim, NULL,PETSC_DECIDE);
1083: PetscOptionsTail();
1084: PetscOptionsEnd();
1085: if (zerobase) shift = 0;
1087: DMCreate(comm, dm);
1088: DMSetType(*dm, DMPLEX);
1089: PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);
1091: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);
1093: /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1094: if (binary) {
1095: parentviewer = viewer;
1096: PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1097: }
1099: if (!rank) {
1100: GmshFile gmsh_, *gmsh = &gmsh_;
1101: char line[PETSC_MAX_PATH_LEN];
1102: PetscBool match;
1104: PetscArrayzero(gmsh,1);
1105: gmsh->viewer = viewer;
1106: gmsh->binary = binary;
1108: /* Read mesh format */
1109: GmshReadSection(gmsh, line);
1110: GmshExpect(gmsh, "$MeshFormat", line);
1111: DMPlexCreateGmsh_ReadMeshFormat(gmsh);
1112: GmshReadEndSection(gmsh, "$EndMeshFormat", line);
1114: /* OPTIONAL Read physical names */
1115: GmshReadSection(gmsh, line);
1116: GmshMatch(gmsh,"$PhysicalNames", line, &match);
1117: if (match) {
1118: DMPlexCreateGmsh_ReadPhysicalNames(gmsh);
1119: GmshReadEndSection(gmsh, "$EndPhysicalNames", line);
1120: /* Initial read for entity section */
1121: GmshReadSection(gmsh, line);
1122: }
1124: /* Read entities */
1125: if (gmsh->fileFormat >= 40) {
1126: GmshExpect(gmsh, "$Entities", line);
1127: switch (gmsh->fileFormat) {
1128: case 41: DMPlexCreateGmsh_ReadEntities_v41(gmsh, &entities); break;
1129: default: DMPlexCreateGmsh_ReadEntities_v40(gmsh, &entities); break;
1130: }
1131: GmshReadEndSection(gmsh, "$EndEntities", line);
1132: /* Initial read for nodes section */
1133: GmshReadSection(gmsh, line);
1134: }
1136: /* Read nodes */
1137: GmshExpect(gmsh, "$Nodes", line);
1138: switch (gmsh->fileFormat) {
1139: case 41: DMPlexCreateGmsh_ReadNodes_v41(gmsh, shift, &numVertices, &coordsIn); break;
1140: case 40: DMPlexCreateGmsh_ReadNodes_v40(gmsh, shift, &numVertices, &coordsIn); break;
1141: default: DMPlexCreateGmsh_ReadNodes_v22(gmsh, shift, &numVertices, &coordsIn); break;
1142: }
1143: GmshReadEndSection(gmsh, "$EndNodes", line);
1145: /* Read elements */
1146: GmshReadSection(gmsh, line);
1147: GmshExpect(gmsh, "$Elements", line);
1148: switch (gmsh->fileFormat) {
1149: case 41: DMPlexCreateGmsh_ReadElements_v41(gmsh, shift, entities, &numCells, &gmsh_elem); break;
1150: case 40: DMPlexCreateGmsh_ReadElements_v40(gmsh, shift, entities, &numCells, &gmsh_elem); break;
1151: default: DMPlexCreateGmsh_ReadElements_v22(gmsh, shift, &numCells, &gmsh_elem); break;
1152: }
1153: GmshReadEndSection(gmsh, "$EndElements", line);
1155: /* OPTIONAL Read periodic section */
1156: if (periodic) {
1157: GmshReadSection(gmsh, line);
1158: GmshMatch(gmsh, "$Periodic", line, &periodic);
1159: }
1160: if (periodic) {
1161: PetscInt pVert, *periodicMapT, *aux;
1163: PetscMalloc1(numVertices, &periodicMapT);
1164: PetscBTCreate(numVertices, &periodicV);
1165: for (i = 0; i < numVertices; i++) periodicMapT[i] = i;
1167: GmshExpect(gmsh, "$Periodic", line);
1168: switch (gmsh->fileFormat) {
1169: case 41: DMPlexCreateGmsh_ReadPeriodic_v41(gmsh, shift, periodicMapT, periodicV); break;
1170: default: DMPlexCreateGmsh_ReadPeriodic_v40(gmsh, shift, periodicMapT, periodicV); break;
1171: }
1172: GmshReadEndSection(gmsh, "$EndPeriodic", line);
1174: /* we may have slaves of slaves */
1175: for (i = 0; i < numVertices; i++) {
1176: while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) {
1177: periodicMapT[i] = periodicMapT[periodicMapT[i]];
1178: }
1179: }
1180: /* periodicMap : from old to new numbering (periodic vertices excluded)
1181: periodicMapI: from new to old numbering */
1182: PetscMalloc1(numVertices, &periodicMap);
1183: PetscMalloc1(numVertices, &periodicMapI);
1184: PetscMalloc1(numVertices, &aux);
1185: for (i = 0, pVert = 0; i < numVertices; i++) {
1186: if (periodicMapT[i] != i) {
1187: pVert++;
1188: } else {
1189: aux[i] = i - pVert;
1190: periodicMapI[i - pVert] = i;
1191: }
1192: }
1193: for (i = 0 ; i < numVertices; i++) {
1194: periodicMap[i] = aux[periodicMapT[i]];
1195: }
1196: PetscFree(periodicMapT);
1197: PetscFree(aux);
1198: /* remove periodic vertices */
1199: numVertices = numVertices - pVert;
1200: }
1202: GmshEntitiesDestroy(&entities);
1203: PetscFree(gmsh->wbuf);
1204: PetscFree(gmsh->sbuf);
1205: }
1207: if (parentviewer) {
1208: PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1209: }
1211: if (!rank) {
1212: PetscBool hybrid = PETSC_FALSE;
1213: PetscInt cellType = -1;
1215: for (trueNumCells = 0, c = 0; c < numCells; ++c) {
1216: if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0; hybrid = PETSC_FALSE; cellType = -1;}
1217: if (gmsh_elem[c].dim < dim) continue;
1218: if (cellType == -1) cellType = gmsh_elem[c].cellType;
1219: /* different cell type indicate an hybrid mesh in PLEX */
1220: if (cellType != gmsh_elem[c].cellType) hybrid = PETSC_TRUE;
1221: /* wedges always indicate an hybrid mesh in PLEX */
1222: if (cellType == 6 || cellType == 13) hybrid = PETSC_TRUE;
1223: if (cellType == 4 || cellType == 11) hasTetra = PETSC_TRUE;
1224: trueNumCells++;
1225: }
1226: /* Renumber cells for hybrid grids */
1227: if (hybrid && enable_hybrid) {
1228: PetscInt hc1 = 0, hc2 = 0, *hybridCells1 = NULL, *hybridCells2 = NULL;
1229: PetscInt cell, tn, *tp;
1230: int n1 = 0,n2 = 0;
1232: PetscMalloc1(trueNumCells, &hybridCells1);
1233: PetscMalloc1(trueNumCells, &hybridCells2);
1234: for (cell = 0, c = 0; c < numCells; ++c) {
1235: if (gmsh_elem[c].dim == dim) {
1236: if (!n1) n1 = gmsh_elem[c].cellType;
1237: else if (!n2 && n1 != gmsh_elem[c].cellType) n2 = gmsh_elem[c].cellType;
1239: if (gmsh_elem[c].cellType == n1) hybridCells1[hc1++] = cell;
1240: else if (gmsh_elem[c].cellType == n2) hybridCells2[hc2++] = cell;
1241: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle more than 2 cell types");
1242: cell++;
1243: }
1244: }
1246: switch (n1) {
1247: case 2: /* triangles */
1248: case 9:
1249: switch (n2) {
1250: case 0: /* single type mesh */
1251: case 3: /* quads */
1252: case 10:
1253: break;
1254: default:
1255: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1256: }
1257: break;
1258: case 3: /* quadrilateral */
1259: case 10:
1260: switch (n2) {
1261: case 0: /* single type mesh */
1262: case 2: /* swap since we list simplices first */
1263: case 9:
1264: tn = hc1;
1265: hc1 = hc2;
1266: hc2 = tn;
1268: tp = hybridCells1;
1269: hybridCells1 = hybridCells2;
1270: hybridCells2 = tp;
1271: break;
1272: default:
1273: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1274: }
1275: break;
1276: case 4: /* tetrahedra */
1277: case 11:
1278: switch (n2) {
1279: case 0: /* single type mesh */
1280: case 6: /* wedges */
1281: case 13:
1282: break;
1283: default:
1284: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1285: }
1286: break;
1287: case 5: /* hexahedra */
1288: case 12:
1289: switch (n2) {
1290: case 0: /* single type mesh */
1291: case 6: /* wedges */
1292: case 13:
1293: break;
1294: default:
1295: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1296: }
1297: break;
1298: case 6: /* wedge */
1299: case 13:
1300: switch (n2) {
1301: case 0: /* single type mesh */
1302: case 4: /* tetrahedra: swap since we list simplices first */
1303: case 11:
1304: case 5: /* hexahedra */
1305: case 12:
1306: tn = hc1;
1307: hc1 = hc2;
1308: hc2 = tn;
1310: tp = hybridCells1;
1311: hybridCells1 = hybridCells2;
1312: hybridCells2 = tp;
1313: break;
1314: default:
1315: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1316: }
1317: break;
1318: default:
1319: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1320: }
1321: PetscMalloc1(trueNumCells, &hybridMap);
1322: for (cell = 0; cell < hc1; cell++) hybridMap[hybridCells1[cell]] = cell;
1323: for (cell = 0; cell < hc2; cell++) hybridMap[hybridCells2[cell]] = cell + hc1;
1324: PetscFree(hybridCells1);
1325: PetscFree(hybridCells2);
1326: }
1328: }
1330: /* Allocate the cell-vertex mesh */
1331: /* We do not want this label automatically computed, instead we compute it here */
1332: DMCreateLabel(*dm, "celltype");
1333: DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
1334: for (cell = 0, c = 0; c < numCells; ++c) {
1335: if (gmsh_elem[c].dim == dim) {
1336: PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1337: DMPolytopeType ctype = DMPolytopeTypeFromGmsh(gmsh_elem[c].cellType);
1338: if (hybridMap && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR;
1339: DMPlexSetConeSize(*dm, ucell, gmsh_elem[c].numNodes);
1340: DMPlexSetCellType(*dm, ucell, ctype);
1341: cell++;
1342: }
1343: }
1344: for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
1345: DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
1346: }
1347: DMSetUp(*dm);
1349: /* Add cell-vertex connections */
1350: for (cell = 0, c = 0; c < numCells; ++c) {
1351: if (gmsh_elem[c].dim == dim) {
1352: PetscInt pcone[8], corner;
1353: PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1354: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1355: const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
1356: pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells;
1357: }
1358: DMPlexReorderCell(*dm, ucell, pcone);
1359: DMPlexSetCone(*dm, ucell, pcone);
1360: cell++;
1361: }
1362: }
1364: MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
1365: DMSetDimension(*dm, dim);
1366: DMPlexSymmetrize(*dm);
1367: DMPlexStratify(*dm);
1368: if (interpolate) {
1369: DM idm;
1371: DMPlexInterpolate(*dm, &idm);
1372: DMDestroy(dm);
1373: *dm = idm;
1374: }
1376: if (usemarker && !interpolate && dim > 1) SETERRQ(comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation");
1377: if (!rank && usemarker) {
1378: PetscInt f, fStart, fEnd;
1380: DMCreateLabel(*dm, "marker");
1381: DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
1382: for (f = fStart; f < fEnd; ++f) {
1383: PetscInt suppSize;
1385: DMPlexGetSupportSize(*dm, f, &suppSize);
1386: if (suppSize == 1) {
1387: PetscInt *cone = NULL, coneSize, p;
1389: DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1390: for (p = 0; p < coneSize; p += 2) {
1391: DMSetLabelValue(*dm, "marker", cone[p], 1);
1392: }
1393: DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1394: }
1395: }
1396: }
1398: if (!rank) {
1399: PetscInt vStart, vEnd;
1401: DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
1402: for (cell = 0, c = 0; c < numCells; ++c) {
1404: /* Create face sets */
1405: if (interpolate && gmsh_elem[c].dim == dim-1) {
1406: const PetscInt *join;
1407: PetscInt joinSize, pcone[8], corner;
1408: /* Find the relevant facet with vertex joins */
1409: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1410: const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
1411: pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart;
1412: }
1413: DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, pcone, &joinSize, &join);
1414: if (joinSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for GMsh element %d (Plex cell %D)", gmsh_elem[c].id, c);
1415: DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
1416: DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
1417: }
1419: /* Create cell sets */
1420: if (gmsh_elem[c].dim == dim) {
1421: if (gmsh_elem[c].numTags > 0) {
1422: DMSetLabelValue(*dm, "Cell Sets", hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].tags[0]);
1423: }
1424: cell++;
1425: }
1427: /* Create vertex sets */
1428: if (gmsh_elem[c].dim == 0) {
1429: if (gmsh_elem[c].numTags > 0) {
1430: const PetscInt cc = gmsh_elem[c].nodes[0] - shift;
1431: const PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart;
1432: DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);
1433: }
1434: }
1435: }
1436: }
1438: /* Create coordinates */
1439: if (embedDim < 0) embedDim = dim;
1440: DMSetCoordinateDim(*dm, embedDim);
1441: DMGetCoordinateSection(*dm, &coordSection);
1442: PetscSectionSetNumFields(coordSection, 1);
1443: PetscSectionSetFieldComponents(coordSection, 0, embedDim);
1444: if (periodicMap) { /* we need to localize coordinates on cells */
1445: PetscSectionSetChart(coordSection, 0, trueNumCells + numVertices);
1446: } else {
1447: PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);
1448: }
1449: for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
1450: PetscSectionSetDof(coordSection, v, embedDim);
1451: PetscSectionSetFieldDof(coordSection, v, 0, embedDim);
1452: }
1453: if (periodicMap) {
1454: PetscBTCreate(trueNumCells, &periodicC);
1455: for (cell = 0, c = 0; c < numCells; ++c) {
1456: if (gmsh_elem[c].dim == dim) {
1457: PetscInt corner, pc = PETSC_FALSE;
1458: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1459: pc = (PetscBool)(pc || PetscBTLookup(periodicV, gmsh_elem[c].nodes[corner] - shift));
1460: }
1461: if (pc) {
1462: PetscInt dof = gmsh_elem[c].numNodes*embedDim;
1463: PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1464: PetscBTSet(periodicC, ucell);
1465: PetscSectionSetDof(coordSection, ucell, dof);
1466: PetscSectionSetFieldDof(coordSection, ucell, 0, dof);
1467: }
1468: cell++;
1469: }
1470: }
1471: }
1472: PetscSectionSetUp(coordSection);
1473: PetscSectionGetStorageSize(coordSection, &coordSize);
1474: VecCreate(PETSC_COMM_SELF, &coordinates);
1475: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1476: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1477: VecSetBlockSize(coordinates, embedDim);
1478: VecSetType(coordinates, VECSTANDARD);
1479: VecGetArray(coordinates, &coords);
1480: if (periodicMap) {
1481: PetscInt off;
1483: for (cell = 0, c = 0; c < numCells; ++c) {
1484: if (gmsh_elem[c].dim == dim) {
1485: PetscInt pcone[8], corner;
1486: PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1487: if (PetscUnlikely(PetscBTLookup(periodicC, ucell))) {
1488: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1489: pcone[corner] = gmsh_elem[c].nodes[corner] - shift;
1490: }
1491: DMPlexReorderCell(*dm, ucell, pcone);
1492: PetscSectionGetOffset(coordSection, ucell, &off);
1493: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner)
1494: for (v = pcone[corner], d = 0; d < embedDim; ++d)
1495: coords[off++] = (PetscReal) coordsIn[v*3+d];
1496: }
1497: cell++;
1498: }
1499: }
1500: for (v = 0; v < numVertices; ++v) {
1501: PetscSectionGetOffset(coordSection, v + trueNumCells, &off);
1502: for (d = 0; d < embedDim; ++d) {
1503: coords[off+d] = (PetscReal) coordsIn[periodicMapI[v]*3+d];
1504: }
1505: }
1506: } else {
1507: for (v = 0; v < numVertices; ++v) {
1508: for (d = 0; d < embedDim; ++d) {
1509: coords[v*embedDim+d] = (PetscReal) coordsIn[v*3+d];
1510: }
1511: }
1512: }
1513: VecRestoreArray(coordinates, &coords);
1514: DMSetCoordinatesLocal(*dm, coordinates);
1516: periodic = periodicMap ? PETSC_TRUE : PETSC_FALSE;
1517: MPI_Bcast(&periodic, 1, MPIU_BOOL, 0, comm);
1518: DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);
1520: PetscFree(coordsIn);
1521: PetscFree(gmsh_elem);
1522: PetscFree(hybridMap);
1523: PetscFree(periodicMap);
1524: PetscFree(periodicMapI);
1525: PetscBTDestroy(&periodicV);
1526: PetscBTDestroy(&periodicC);
1527: VecDestroy(&coordinates);
1529: PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
1530: return(0);
1531: }