Actual source code: plexgmsh.c
petsc-3.12.5 2020-03-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: }
943: /*@C
944: DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file
946: + comm - The MPI communicator
947: . filename - Name of the Gmsh file
948: - interpolate - Create faces and edges in the mesh
950: Output Parameter:
951: . dm - The DM object representing the mesh
953: Level: beginner
955: .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
956: @*/
957: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
958: {
959: PetscViewer viewer;
960: PetscMPIInt rank;
961: int fileType;
962: PetscViewerType vtype;
963: PetscErrorCode ierr;
966: MPI_Comm_rank(comm, &rank);
968: /* Determine Gmsh file type (ASCII or binary) from file header */
969: if (!rank) {
970: GmshFile gmsh_, *gmsh = &gmsh_;
971: char line[PETSC_MAX_PATH_LEN];
972: int snum;
973: float version;
975: PetscArrayzero(gmsh,1);
976: PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer);
977: PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII);
978: PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ);
979: PetscViewerFileSetName(gmsh->viewer, filename);
980: /* Read only the first two lines of the Gmsh file */
981: GmshReadSection(gmsh, line);
982: GmshExpect(gmsh, "$MeshFormat", line);
983: GmshReadString(gmsh, line, 2);
984: snum = sscanf(line, "%f %d", &version, &fileType);
985: if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
986: 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);
987: if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f not supported", (double)version);
988: 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);
989: PetscViewerDestroy(&gmsh->viewer);
990: }
991: MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);
992: vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;
994: /* Create appropriate viewer and build plex */
995: PetscViewerCreate(comm, &viewer);
996: PetscViewerSetType(viewer, vtype);
997: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
998: PetscViewerFileSetName(viewer, filename);
999: DMPlexCreateGmsh(comm, viewer, interpolate, dm);
1000: PetscViewerDestroy(&viewer);
1001: return(0);
1002: }
1004: /*@
1005: DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer
1007: Collective
1009: Input Parameters:
1010: + comm - The MPI communicator
1011: . viewer - The Viewer associated with a Gmsh file
1012: - interpolate - Create faces and edges in the mesh
1014: Output Parameter:
1015: . dm - The DM object representing the mesh
1017: Note: http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
1019: Level: beginner
1021: .seealso: DMPLEX, DMCreate()
1022: @*/
1023: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1024: {
1025: PetscViewer parentviewer = NULL;
1026: double *coordsIn = NULL;
1027: GmshEntities *entities = NULL;
1028: GmshElement *gmsh_elem = NULL;
1029: PetscSection coordSection;
1030: Vec coordinates;
1031: PetscBT periodicV = NULL, periodicC = NULL;
1032: PetscScalar *coords;
1033: PetscInt dim = 0, embedDim = -1, coordSize, c, v, d, cell, *periodicMap = NULL, *periodicMapI = NULL, *hybridMap = NULL, cMax = PETSC_DETERMINE;
1034: PetscInt numVertices = 0, numCells = 0, trueNumCells = 0;
1035: int i, shift = 1;
1036: PetscMPIInt rank;
1037: PetscBool binary, zerobase = PETSC_FALSE, usemarker = PETSC_FALSE;
1038: PetscBool enable_hybrid = interpolate, periodic = PETSC_TRUE;
1042: MPI_Comm_rank(comm, &rank);
1043: PetscObjectOptionsBegin((PetscObject)viewer);
1044: PetscOptionsHead(PetscOptionsObject,"DMPlex Gmsh options");
1045: PetscOptionsBool("-dm_plex_gmsh_hybrid", "Generate hybrid cell bounds", "DMPlexCreateGmsh", enable_hybrid, &enable_hybrid, NULL);
1046: PetscOptionsBool("-dm_plex_gmsh_periodic","Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL);
1047: PetscOptionsBool("-dm_plex_gmsh_use_marker", "Generate marker label", "DMPlexCreateGmsh", usemarker, &usemarker, NULL);
1048: PetscOptionsBool("-dm_plex_gmsh_zero_base", "Read Gmsh file with zero base indices", "DMPlexCreateGmsh", zerobase, &zerobase, NULL);
1049: PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", embedDim, &embedDim, NULL,PETSC_DECIDE);
1050: PetscOptionsTail();
1051: PetscOptionsEnd();
1052: if (zerobase) shift = 0;
1054: DMCreate(comm, dm);
1055: DMSetType(*dm, DMPLEX);
1056: PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);
1058: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);
1060: /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1061: if (binary) {
1062: parentviewer = viewer;
1063: PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1064: }
1066: if (!rank) {
1067: GmshFile gmsh_, *gmsh = &gmsh_;
1068: char line[PETSC_MAX_PATH_LEN];
1069: PetscBool match;
1071: PetscArrayzero(gmsh,1);
1072: gmsh->viewer = viewer;
1073: gmsh->binary = binary;
1075: /* Read mesh format */
1076: GmshReadSection(gmsh, line);
1077: GmshExpect(gmsh, "$MeshFormat", line);
1078: DMPlexCreateGmsh_ReadMeshFormat(gmsh);
1079: GmshReadEndSection(gmsh, "$EndMeshFormat", line);
1081: /* OPTIONAL Read physical names */
1082: GmshReadSection(gmsh, line);
1083: GmshMatch(gmsh,"$PhysicalNames", line, &match);
1084: if (match) {
1085: DMPlexCreateGmsh_ReadPhysicalNames(gmsh);
1086: GmshReadEndSection(gmsh, "$EndPhysicalNames", line);
1087: /* Initial read for entity section */
1088: GmshReadSection(gmsh, line);
1089: }
1091: /* Read entities */
1092: if (gmsh->fileFormat >= 40) {
1093: GmshExpect(gmsh, "$Entities", line);
1094: switch (gmsh->fileFormat) {
1095: case 41: DMPlexCreateGmsh_ReadEntities_v41(gmsh, &entities); break;
1096: default: DMPlexCreateGmsh_ReadEntities_v40(gmsh, &entities); break;
1097: }
1098: GmshReadEndSection(gmsh, "$EndEntities", line);
1099: /* Initial read for nodes section */
1100: GmshReadSection(gmsh, line);
1101: }
1103: /* Read nodes */
1104: GmshExpect(gmsh, "$Nodes", line);
1105: switch (gmsh->fileFormat) {
1106: case 41: DMPlexCreateGmsh_ReadNodes_v41(gmsh, shift, &numVertices, &coordsIn); break;
1107: case 40: DMPlexCreateGmsh_ReadNodes_v40(gmsh, shift, &numVertices, &coordsIn); break;
1108: default: DMPlexCreateGmsh_ReadNodes_v22(gmsh, shift, &numVertices, &coordsIn); break;
1109: }
1110: GmshReadEndSection(gmsh, "$EndNodes", line);
1112: /* Read elements */
1113: GmshReadSection(gmsh, line);
1114: GmshExpect(gmsh, "$Elements", line);
1115: switch (gmsh->fileFormat) {
1116: case 41: DMPlexCreateGmsh_ReadElements_v41(gmsh, shift, entities, &numCells, &gmsh_elem); break;
1117: case 40: DMPlexCreateGmsh_ReadElements_v40(gmsh, shift, entities, &numCells, &gmsh_elem); break;
1118: default: DMPlexCreateGmsh_ReadElements_v22(gmsh, shift, &numCells, &gmsh_elem); break;
1119: }
1120: GmshReadEndSection(gmsh, "$EndElements", line);
1122: /* OPTIONAL Read periodic section */
1123: if (periodic) {
1124: GmshReadSection(gmsh, line);
1125: GmshMatch(gmsh, "$Periodic", line, &periodic);
1126: }
1127: if (periodic) {
1128: PetscInt pVert, *periodicMapT, *aux;
1130: PetscMalloc1(numVertices, &periodicMapT);
1131: PetscBTCreate(numVertices, &periodicV);
1132: for (i = 0; i < numVertices; i++) periodicMapT[i] = i;
1134: GmshExpect(gmsh, "$Periodic", line);
1135: switch (gmsh->fileFormat) {
1136: case 41: DMPlexCreateGmsh_ReadPeriodic_v41(gmsh, shift, periodicMapT, periodicV); break;
1137: default: DMPlexCreateGmsh_ReadPeriodic_v40(gmsh, shift, periodicMapT, periodicV); break;
1138: }
1139: GmshReadEndSection(gmsh, "$EndPeriodic", line);
1141: /* we may have slaves of slaves */
1142: for (i = 0; i < numVertices; i++) {
1143: while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) {
1144: periodicMapT[i] = periodicMapT[periodicMapT[i]];
1145: }
1146: }
1147: /* periodicMap : from old to new numbering (periodic vertices excluded)
1148: periodicMapI: from new to old numbering */
1149: PetscMalloc1(numVertices, &periodicMap);
1150: PetscMalloc1(numVertices, &periodicMapI);
1151: PetscMalloc1(numVertices, &aux);
1152: for (i = 0, pVert = 0; i < numVertices; i++) {
1153: if (periodicMapT[i] != i) {
1154: pVert++;
1155: } else {
1156: aux[i] = i - pVert;
1157: periodicMapI[i - pVert] = i;
1158: }
1159: }
1160: for (i = 0 ; i < numVertices; i++) {
1161: periodicMap[i] = aux[periodicMapT[i]];
1162: }
1163: PetscFree(periodicMapT);
1164: PetscFree(aux);
1165: /* remove periodic vertices */
1166: numVertices = numVertices - pVert;
1167: }
1169: GmshEntitiesDestroy(&entities);
1170: PetscFree(gmsh->wbuf);
1171: PetscFree(gmsh->sbuf);
1172: }
1174: if (parentviewer) {
1175: PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1176: }
1178: if (!rank) {
1179: PetscBool hybrid = PETSC_FALSE;
1180: PetscInt cellType = -1;
1182: for (trueNumCells = 0, c = 0; c < numCells; ++c) {
1183: if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0; hybrid = PETSC_FALSE; cellType = -1;}
1184: if (gmsh_elem[c].dim < dim) continue;
1185: if (cellType == -1) cellType = gmsh_elem[c].cellType;
1186: /* different cell type indicate an hybrid mesh in PLEX */
1187: if (cellType != gmsh_elem[c].cellType) hybrid = PETSC_TRUE;
1188: /* wedges always indicate an hybrid mesh in PLEX */
1189: if (cellType == 6 || cellType == 13) hybrid = PETSC_TRUE;
1190: trueNumCells++;
1191: }
1192: /* Renumber cells for hybrid grids */
1193: if (hybrid && enable_hybrid) {
1194: PetscInt hc1 = 0, hc2 = 0, *hybridCells1 = NULL, *hybridCells2 = NULL;
1195: PetscInt cell, tn, *tp;
1196: int n1 = 0,n2 = 0;
1198: PetscMalloc1(trueNumCells, &hybridCells1);
1199: PetscMalloc1(trueNumCells, &hybridCells2);
1200: for (cell = 0, c = 0; c < numCells; ++c) {
1201: if (gmsh_elem[c].dim == dim) {
1202: if (!n1) n1 = gmsh_elem[c].cellType;
1203: else if (!n2 && n1 != gmsh_elem[c].cellType) n2 = gmsh_elem[c].cellType;
1205: if (gmsh_elem[c].cellType == n1) hybridCells1[hc1++] = cell;
1206: else if (gmsh_elem[c].cellType == n2) hybridCells2[hc2++] = cell;
1207: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle more than 2 cell types");
1208: cell++;
1209: }
1210: }
1212: switch (n1) {
1213: case 2: /* triangles */
1214: case 9:
1215: switch (n2) {
1216: case 0: /* single type mesh */
1217: case 3: /* quads */
1218: case 10:
1219: break;
1220: default:
1221: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1222: }
1223: break;
1224: case 3: /* quadrilateral */
1225: case 10:
1226: switch (n2) {
1227: case 0: /* single type mesh */
1228: case 2: /* swap since we list simplices first */
1229: case 9:
1230: tn = hc1;
1231: hc1 = hc2;
1232: hc2 = tn;
1234: tp = hybridCells1;
1235: hybridCells1 = hybridCells2;
1236: hybridCells2 = tp;
1237: break;
1238: default:
1239: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1240: }
1241: break;
1242: case 4: /* tetrahedra */
1243: case 11:
1244: switch (n2) {
1245: case 0: /* single type mesh */
1246: case 6: /* wedges */
1247: case 13:
1248: break;
1249: default:
1250: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1251: }
1252: break;
1253: case 5: /* hexahedra */
1254: case 12:
1255: switch (n2) {
1256: case 0: /* single type mesh */
1257: case 6: /* wedges */
1258: case 13:
1259: break;
1260: default:
1261: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1262: }
1263: break;
1264: case 6: /* wedge */
1265: case 13:
1266: switch (n2) {
1267: case 0: /* single type mesh */
1268: case 4: /* tetrahedra: swap since we list simplices first */
1269: case 11:
1270: case 5: /* hexahedra */
1271: case 12:
1272: tn = hc1;
1273: hc1 = hc2;
1274: hc2 = tn;
1276: tp = hybridCells1;
1277: hybridCells1 = hybridCells2;
1278: hybridCells2 = tp;
1279: break;
1280: default:
1281: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1282: }
1283: break;
1284: default:
1285: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1286: }
1287: cMax = hc1;
1288: PetscMalloc1(trueNumCells, &hybridMap);
1289: for (cell = 0; cell < hc1; cell++) hybridMap[hybridCells1[cell]] = cell;
1290: for (cell = 0; cell < hc2; cell++) hybridMap[hybridCells2[cell]] = cell + hc1;
1291: PetscFree(hybridCells1);
1292: PetscFree(hybridCells2);
1293: }
1295: }
1297: /* Allocate the cell-vertex mesh */
1298: DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
1299: for (cell = 0, c = 0; c < numCells; ++c) {
1300: if (gmsh_elem[c].dim == dim) {
1301: DMPlexSetConeSize(*dm, hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].numNodes);
1302: cell++;
1303: }
1304: }
1305: DMSetUp(*dm);
1306: /* Add cell-vertex connections */
1307: for (cell = 0, c = 0; c < numCells; ++c) {
1308: if (gmsh_elem[c].dim == dim) {
1309: PetscInt pcone[8], corner;
1310: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1311: const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
1312: pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells;
1313: }
1314: if (dim == 3) {
1315: /* Tetrahedra are inverted */
1316: if (gmsh_elem[c].cellType == 4 || gmsh_elem[c].cellType == 11) {
1317: PetscInt tmp = pcone[0];
1318: pcone[0] = pcone[1];
1319: pcone[1] = tmp;
1320: }
1321: /* Hexahedra are inverted */
1322: if (gmsh_elem[c].cellType == 5 || gmsh_elem[c].cellType == 12) {
1323: PetscInt tmp = pcone[1];
1324: pcone[1] = pcone[3];
1325: pcone[3] = tmp;
1326: }
1327: /* Prisms are inverted */
1328: if (gmsh_elem[c].cellType == 6 || gmsh_elem[c].cellType == 13) {
1329: PetscInt tmp;
1331: tmp = pcone[1];
1332: pcone[1] = pcone[2];
1333: pcone[2] = tmp;
1334: tmp = pcone[4];
1335: pcone[4] = pcone[5];
1336: pcone[5] = tmp;
1337: }
1338: } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */
1339: /* quads are input to PLEX as prisms */
1340: if (gmsh_elem[c].cellType == 3 || gmsh_elem[c].cellType == 10) {
1341: PetscInt tmp = pcone[2];
1342: pcone[2] = pcone[3];
1343: pcone[3] = tmp;
1344: }
1345: }
1346: DMPlexSetCone(*dm, hybridMap ? hybridMap[cell] : cell, pcone);
1347: cell++;
1348: }
1349: }
1350: MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
1351: DMSetDimension(*dm, dim);
1352: DMPlexSetHybridBounds(*dm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
1353: DMPlexSymmetrize(*dm);
1354: DMPlexStratify(*dm);
1355: if (interpolate) {
1356: DM idm;
1358: DMPlexInterpolate(*dm, &idm);
1359: DMDestroy(dm);
1360: *dm = idm;
1361: }
1363: if (usemarker && !interpolate && dim > 1) SETERRQ(comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation");
1364: if (!rank && usemarker) {
1365: PetscInt f, fStart, fEnd;
1367: DMCreateLabel(*dm, "marker");
1368: DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
1369: for (f = fStart; f < fEnd; ++f) {
1370: PetscInt suppSize;
1372: DMPlexGetSupportSize(*dm, f, &suppSize);
1373: if (suppSize == 1) {
1374: PetscInt *cone = NULL, coneSize, p;
1376: DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1377: for (p = 0; p < coneSize; p += 2) {
1378: DMSetLabelValue(*dm, "marker", cone[p], 1);
1379: }
1380: DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1381: }
1382: }
1383: }
1385: if (!rank) {
1386: PetscInt vStart, vEnd;
1388: DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
1389: for (cell = 0, c = 0; c < numCells; ++c) {
1391: /* Create face sets */
1392: if (interpolate && gmsh_elem[c].dim == dim-1) {
1393: const PetscInt *join;
1394: PetscInt joinSize, pcone[8], corner;
1395: /* Find the relevant facet with vertex joins */
1396: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1397: const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
1398: pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart;
1399: }
1400: DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, pcone, &joinSize, &join);
1401: 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);
1402: DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
1403: DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
1404: }
1406: /* Create cell sets */
1407: if (gmsh_elem[c].dim == dim) {
1408: if (gmsh_elem[c].numTags > 0) {
1409: DMSetLabelValue(*dm, "Cell Sets", hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].tags[0]);
1410: }
1411: cell++;
1412: }
1414: /* Create vertex sets */
1415: if (gmsh_elem[c].dim == 0) {
1416: if (gmsh_elem[c].numTags > 0) {
1417: const PetscInt cc = gmsh_elem[c].nodes[0] - shift;
1418: const PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart;
1419: DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);
1420: }
1421: }
1422: }
1423: }
1425: /* Create coordinates */
1426: if (embedDim < 0) embedDim = dim;
1427: DMSetCoordinateDim(*dm, embedDim);
1428: DMGetCoordinateSection(*dm, &coordSection);
1429: PetscSectionSetNumFields(coordSection, 1);
1430: PetscSectionSetFieldComponents(coordSection, 0, embedDim);
1431: if (periodicMap) { /* we need to localize coordinates on cells */
1432: PetscSectionSetChart(coordSection, 0, trueNumCells + numVertices);
1433: } else {
1434: PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);
1435: }
1436: for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
1437: PetscSectionSetDof(coordSection, v, embedDim);
1438: PetscSectionSetFieldDof(coordSection, v, 0, embedDim);
1439: }
1440: if (periodicMap) {
1441: PetscBTCreate(trueNumCells, &periodicC);
1442: for (cell = 0, c = 0; c < numCells; ++c) {
1443: if (gmsh_elem[c].dim == dim) {
1444: PetscInt corner;
1445: PetscBool pc = PETSC_FALSE;
1446: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1447: pc = (PetscBool)(pc || PetscBTLookup(periodicV, gmsh_elem[c].nodes[corner] - shift));
1448: }
1449: if (pc) {
1450: PetscInt dof = gmsh_elem[c].numNodes*embedDim;
1451: PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1452: PetscBTSet(periodicC, ucell);
1453: PetscSectionSetDof(coordSection, ucell, dof);
1454: PetscSectionSetFieldDof(coordSection, ucell, 0, dof);
1455: }
1456: cell++;
1457: }
1458: }
1459: }
1460: PetscSectionSetUp(coordSection);
1461: PetscSectionGetStorageSize(coordSection, &coordSize);
1462: VecCreate(PETSC_COMM_SELF, &coordinates);
1463: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1464: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1465: VecSetBlockSize(coordinates, embedDim);
1466: VecSetType(coordinates, VECSTANDARD);
1467: VecGetArray(coordinates, &coords);
1468: if (periodicMap) {
1469: PetscInt off;
1471: for (cell = 0, c = 0; c < numCells; ++c) {
1472: PetscInt pcone[8], corner;
1473: if (gmsh_elem[c].dim == dim) {
1474: PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1475: if (PetscUnlikely(PetscBTLookup(periodicC, ucell))) {
1476: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1477: pcone[corner] = gmsh_elem[c].nodes[corner] - shift;
1478: }
1479: if (dim == 3) {
1480: /* Tetrahedra are inverted */
1481: if (gmsh_elem[c].cellType == 4 || gmsh_elem[c].cellType == 11) {
1482: PetscInt tmp = pcone[0];
1483: pcone[0] = pcone[1];
1484: pcone[1] = tmp;
1485: }
1486: /* Hexahedra are inverted */
1487: if (gmsh_elem[c].cellType == 5 || gmsh_elem[c].cellType == 12) {
1488: PetscInt tmp = pcone[1];
1489: pcone[1] = pcone[3];
1490: pcone[3] = tmp;
1491: }
1492: /* Prisms are inverted */
1493: if (gmsh_elem[c].cellType == 6 || gmsh_elem[c].cellType == 13) {
1494: PetscInt tmp;
1496: tmp = pcone[1];
1497: pcone[1] = pcone[2];
1498: pcone[2] = tmp;
1499: tmp = pcone[4];
1500: pcone[4] = pcone[5];
1501: pcone[5] = tmp;
1502: }
1503: } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */
1504: /* quads are input to PLEX as prisms */
1505: if (gmsh_elem[c].cellType == 3 || gmsh_elem[c].cellType == 10) {
1506: PetscInt tmp = pcone[2];
1507: pcone[2] = pcone[3];
1508: pcone[3] = tmp;
1509: }
1510: }
1511: PetscSectionGetOffset(coordSection, ucell, &off);
1512: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1513: v = pcone[corner];
1514: for (d = 0; d < embedDim; ++d) {
1515: coords[off++] = (PetscReal) coordsIn[v*3+d];
1516: }
1517: }
1518: }
1519: cell++;
1520: }
1521: }
1522: for (v = 0; v < numVertices; ++v) {
1523: PetscSectionGetOffset(coordSection, v + trueNumCells, &off);
1524: for (d = 0; d < embedDim; ++d) {
1525: coords[off+d] = (PetscReal) coordsIn[periodicMapI[v]*3+d];
1526: }
1527: }
1528: } else {
1529: for (v = 0; v < numVertices; ++v) {
1530: for (d = 0; d < embedDim; ++d) {
1531: coords[v*embedDim+d] = (PetscReal) coordsIn[v*3+d];
1532: }
1533: }
1534: }
1535: VecRestoreArray(coordinates, &coords);
1536: DMSetCoordinatesLocal(*dm, coordinates);
1538: periodic = periodicMap ? PETSC_TRUE : PETSC_FALSE;
1539: MPI_Bcast(&periodic, 1, MPIU_BOOL, 0, comm);
1540: DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);
1542: PetscFree(coordsIn);
1543: PetscFree(gmsh_elem);
1544: PetscFree(hybridMap);
1545: PetscFree(periodicMap);
1546: PetscFree(periodicMapI);
1547: PetscBTDestroy(&periodicV);
1548: PetscBTDestroy(&periodicC);
1549: VecDestroy(&coordinates);
1551: PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
1552: return(0);
1553: }