Actual source code: plexgmsh.c
petsc-3.11.4 2019-09-28
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 13: /* 18-node 2nd wedge */
330: dim = 3;
331: numNodes = 6;
332: numNodesIgnore = 12;
333: break;
334: case 15: /* 1-node vertex */
335: dim = 0;
336: numNodes = 1;
337: break;
338: case 7: /* 5-node pyramid */
339: case 10: /* 9-node 2nd order quadrangle */
340: case 11: /* 10-node 2nd order tetrahedron */
341: case 12: /* 27-node 2nd order hexhedron */
342: case 14: /* 14-node 2nd order pyramid */
343: default:
344: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
345: }
346: if (binary) {
347: const int nint = 1 + numTags + numNodes + numNodesIgnore;
348: /* Loop over element blocks */
349: for (i = 0; i < numElem; ++i, ++c) {
350: PetscViewerRead(viewer, ibuf, nint, NULL, PETSC_ENUM);
351: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, nint);}
352: elements[c].dim = dim;
353: elements[c].numNodes = numNodes;
354: elements[c].numTags = numTags;
355: elements[c].id = ibuf[0];
356: elements[c].cellType = cellType;
357: for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p];
358: for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
359: }
360: } else {
361: const int nint = numTags + numNodes + numNodesIgnore;
362: elements[c].dim = dim;
363: elements[c].numNodes = numNodes;
364: elements[c].numTags = numTags;
365: elements[c].cellType = cellType;
366: PetscViewerRead(viewer, ibuf, nint, NULL, PETSC_ENUM);
367: for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[p];
368: for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[numTags + p];
369: c++;
370: }
371: }
372: return(0);
373: }
375: /*
376: $Entities
377: numPoints(unsigned long) numCurves(unsigned long)
378: numSurfaces(unsigned long) numVolumes(unsigned long)
379: // points
380: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
381: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
382: numPhysicals(unsigned long) phyisicalTag[...](int)
383: ...
384: // curves
385: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
386: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
387: numPhysicals(unsigned long) physicalTag[...](int)
388: numBREPVert(unsigned long) tagBREPVert[...](int)
389: ...
390: // surfaces
391: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
392: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
393: numPhysicals(unsigned long) physicalTag[...](int)
394: numBREPCurve(unsigned long) tagBREPCurve[...](int)
395: ...
396: // volumes
397: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
398: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
399: numPhysicals(unsigned long) physicalTag[...](int)
400: numBREPSurfaces(unsigned long) tagBREPSurfaces[...](int)
401: ...
402: $EndEntities
403: */
404: static PetscErrorCode DMPlexCreateGmsh_ReadEntities_v40(GmshFile *gmsh, GmshEntities **entities)
405: {
406: PetscViewer viewer = gmsh->viewer;
407: PetscBool byteSwap = gmsh->byteSwap;
408: long index, num, lbuf[4];
409: int dim, eid, numTags, *ibuf, t;
410: PetscInt count[4], i;
411: GmshEntity *entity = NULL;
415: PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG);
416: if (byteSwap) {PetscByteSwap(lbuf, PETSC_LONG, 4);}
417: for (i = 0; i < 4; ++i) count[i] = lbuf[i];
418: GmshEntitiesCreate(count, entities);
419: for (dim = 0; dim < 4; ++dim) {
420: for (index = 0; index < count[dim]; ++index) {
421: PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM);
422: if (byteSwap) {PetscByteSwap(&eid, PETSC_ENUM, 1);}
423: GmshEntitiesAdd(*entities, (PetscInt)index, dim, eid, &entity);
424: PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE);
425: if (byteSwap) {PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6);}
426: PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
427: if (byteSwap) {PetscByteSwap(&num, PETSC_LONG, 1);}
428: GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
429: PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
430: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, num);}
431: entity->numTags = numTags = (int) PetscMin(num, 4);
432: for (t = 0; t < numTags; ++t) entity->tags[t] = ibuf[t];
433: if (dim == 0) continue;
434: PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
435: if (byteSwap) {PetscByteSwap(&num, PETSC_LONG, 1);}
436: GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
437: PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
438: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, num);}
439: }
440: }
441: return(0);
442: }
444: /*
445: $Nodes
446: numEntityBlocks(unsigned long) numNodes(unsigned long)
447: tagEntity(int) dimEntity(int) typeNode(int) numNodes(unsigned long)
448: tag(int) x(double) y(double) z(double)
449: ...
450: ...
451: $EndNodes
452: */
453: static PetscErrorCode DMPlexCreateGmsh_ReadNodes_v40(GmshFile *gmsh, int shift, PetscInt *numVertices, double **gmsh_nodes)
454: {
455: PetscViewer viewer = gmsh->viewer;
456: PetscBool byteSwap = gmsh->byteSwap;
457: long block, node, v, numEntityBlocks, numTotalNodes, numNodes;
458: int info[3], nid;
459: double *coordinates;
460: char *cbuf;
464: PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
465: if (byteSwap) {PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);}
466: PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG);
467: if (byteSwap) {PetscByteSwap(&numTotalNodes, PETSC_LONG, 1);}
468: PetscMalloc1(numTotalNodes*3, &coordinates);
469: *numVertices = numTotalNodes;
470: *gmsh_nodes = coordinates;
471: for (v = 0, block = 0; block < numEntityBlocks; ++block) {
472: PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
473: PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG);
474: if (byteSwap) {PetscByteSwap(&numNodes, PETSC_LONG, 1);}
475: if (gmsh->binary) {
476: int nbytes = sizeof(int) + 3*sizeof(double);
477: GmshBufferGet(gmsh, numNodes, nbytes, &cbuf);
478: PetscViewerRead(viewer, cbuf, numNodes*nbytes, NULL, PETSC_CHAR);
479: for (node = 0; node < numNodes; ++node, ++v) {
480: char *cnid = cbuf + node*nbytes, *cxyz = cnid + sizeof(int);
481: double *xyz = coordinates + v*3;
482: #if !defined(PETSC_WORDS_BIGENDIAN)
483: PetscByteSwap(cnid, PETSC_ENUM, 1);
484: PetscByteSwap(cxyz, PETSC_DOUBLE, 3);
485: #endif
486: PetscMemcpy(&nid, cnid, sizeof(int));
487: PetscMemcpy(xyz, cxyz, 3*sizeof(double));
488: if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
489: if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
490: if (nid != v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %d", nid, v+shift);
491: }
492: } else {
493: for (node = 0; node < numNodes; ++node, ++v) {
494: double *xyz = coordinates + v*3;
495: PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
496: if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
497: if (nid != v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %d", nid, v+shift);
498: PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
499: if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
500: }
501: }
502: }
503: return(0);
504: }
506: /*
507: $Elements
508: numEntityBlocks(unsigned long) numElements(unsigned long)
509: tagEntity(int) dimEntity(int) typeEle(int) numElements(unsigned long)
510: tag(int) numVert[...](int)
511: ...
512: ...
513: $EndElements
514: */
515: static PetscErrorCode DMPlexCreateGmsh_ReadElements_v40(GmshFile *gmsh, PETSC_UNUSED int shift, GmshEntities *entities, PetscInt *numCells, GmshElement **gmsh_elems)
516: {
517: PetscViewer viewer = gmsh->viewer;
518: PetscBool byteSwap = gmsh->byteSwap;
519: long c, block, numEntityBlocks, numTotalElements, elem, numElements;
520: int p, info[3], *ibuf = NULL;
521: int eid, dim, numTags, *tags, cellType, numNodes;
522: GmshEntity *entity = NULL;
523: GmshElement *elements;
527: PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
528: if (byteSwap) {PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);}
529: PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG);
530: if (byteSwap) {PetscByteSwap(&numTotalElements, PETSC_LONG, 1);}
531: PetscCalloc1(numTotalElements, &elements);
532: *numCells = numTotalElements;
533: *gmsh_elems = elements;
534: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
535: PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
536: if (byteSwap) {PetscByteSwap(info, PETSC_ENUM, 3);}
537: eid = info[0]; dim = info[1]; cellType = info[2];
538: GmshEntitiesGet(entities, dim, eid, &entity);
539: numTags = entity->numTags;
540: tags = entity->tags;
541: switch (cellType) {
542: case 1: /* 2-node line */
543: numNodes = 2;
544: break;
545: case 2: /* 3-node triangle */
546: numNodes = 3;
547: break;
548: case 3: /* 4-node quadrangle */
549: numNodes = 4;
550: break;
551: case 4: /* 4-node tetrahedron */
552: numNodes = 4;
553: break;
554: case 5: /* 8-node hexahedron */
555: numNodes = 8;
556: break;
557: case 6: /* 6-node wedge */
558: numNodes = 6;
559: break;
560: case 15: /* 1-node vertex */
561: numNodes = 1;
562: break;
563: default:
564: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
565: }
566: PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG);
567: if (byteSwap) {PetscByteSwap(&numElements, PETSC_LONG, 1);}
568: GmshBufferGet(gmsh, (1+numNodes)*numElements, sizeof(int), &ibuf);
569: PetscViewerRead(viewer, ibuf, (1+numNodes)*numElements, NULL, PETSC_ENUM);
570: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, (1+numNodes)*numElements);}
571: for (elem = 0; elem < numElements; ++elem, ++c) {
572: int *id = ibuf + elem*(1+numNodes), *nodes = id + 1;
573: GmshElement *element = elements + c;
574: element->dim = dim;
575: element->cellType = cellType;
576: element->numNodes = numNodes;
577: element->numTags = numTags;
578: element->id = id[0];
579: for (p = 0; p < numNodes; p++) element->nodes[p] = nodes[p];
580: for (p = 0; p < numTags; p++) element->tags[p] = tags[p];
581: }
582: }
583: return(0);
584: }
586: static PetscErrorCode DMPlexCreateGmsh_ReadPeriodic_v40(GmshFile *gmsh, int shift, PetscInt slaveMap[], PetscBT bt)
587: {
588: PetscViewer viewer = gmsh->viewer;
589: int fileFormat = gmsh->fileFormat;
590: PetscBool binary = gmsh->binary;
591: PetscBool byteSwap = gmsh->byteSwap;
592: int numPeriodic, snum, i;
593: char line[PETSC_MAX_PATH_LEN];
597: if (fileFormat == 22 || !binary) {
598: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
599: snum = sscanf(line, "%d", &numPeriodic);
600: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
601: } else {
602: PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM);
603: if (byteSwap) {PetscByteSwap(&numPeriodic, PETSC_ENUM, 1);}
604: }
605: for (i = 0; i < numPeriodic; i++) {
606: int ibuf[3], slaveDim = -1, slaveTag = -1, masterTag = -1, slaveNode, masterNode;
607: long j, nNodes;
608: double affine[16];
610: if (fileFormat == 22 || !binary) {
611: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
612: snum = sscanf(line, "%d %d %d", &slaveDim, &slaveTag, &masterTag);
613: if (snum != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
614: } else {
615: PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
616: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 3);}
617: slaveDim = ibuf[0]; slaveTag = ibuf[1]; masterTag = ibuf[2];
618: }
619: (void)slaveDim; (void)slaveTag; (void)masterTag; /* unused */
621: if (fileFormat == 22 || !binary) {
622: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
623: snum = sscanf(line, "%ld", &nNodes);
624: if (snum != 1) { /* discard tranformation and try again */
625: PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
626: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
627: snum = sscanf(line, "%ld", &nNodes);
628: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
629: }
630: } else {
631: PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
632: if (byteSwap) {PetscByteSwap(&nNodes, PETSC_LONG, 1);}
633: if (nNodes == -1) { /* discard tranformation and try again */
634: PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE);
635: PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
636: if (byteSwap) {PetscByteSwap(&nNodes, PETSC_LONG, 1);}
637: }
638: }
640: for (j = 0; j < nNodes; j++) {
641: if (fileFormat == 22 || !binary) {
642: PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
643: snum = sscanf(line, "%d %d", &slaveNode, &masterNode);
644: if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
645: } else {
646: PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM);
647: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 2);}
648: slaveNode = ibuf[0]; masterNode = ibuf[1];
649: }
650: slaveMap[slaveNode - shift] = masterNode - shift;
651: PetscBTSet(bt, slaveNode - shift);
652: PetscBTSet(bt, masterNode - shift);
653: }
654: }
655: return(0);
656: }
658: /*
659: $Entities
660: numPoints(size_t) numCurves(size_t)
661: numSurfaces(size_t) numVolumes(size_t)
662: pointTag(int) X(double) Y(double) Z(double)
663: numPhysicalTags(size_t) physicalTag(int) ...
664: ...
665: curveTag(int) minX(double) minY(double) minZ(double)
666: maxX(double) maxY(double) maxZ(double)
667: numPhysicalTags(size_t) physicalTag(int) ...
668: numBoundingPoints(size_t) pointTag(int) ...
669: ...
670: surfaceTag(int) minX(double) minY(double) minZ(double)
671: maxX(double) maxY(double) maxZ(double)
672: numPhysicalTags(size_t) physicalTag(int) ...
673: numBoundingCurves(size_t) curveTag(int) ...
674: ...
675: volumeTag(int) minX(double) minY(double) minZ(double)
676: maxX(double) maxY(double) maxZ(double)
677: numPhysicalTags(size_t) physicalTag(int) ...
678: numBoundngSurfaces(size_t) surfaceTag(int) ...
679: ...
680: $EndEntities
681: */
682: static PetscErrorCode DMPlexCreateGmsh_ReadEntities_v41(GmshFile *gmsh, GmshEntities **entities)
683: {
684: PetscInt count[4], index, numTags, i;
685: int dim, eid, *tags = NULL;
686: GmshEntity *entity = NULL;
690: GmshReadSize(gmsh, count, 4);
691: GmshEntitiesCreate(count, entities);
692: for (dim = 0; dim < 4; ++dim) {
693: for (index = 0; index < count[dim]; ++index) {
694: GmshReadInt(gmsh, &eid, 1);
695: GmshEntitiesAdd(*entities, (PetscInt)index, dim, eid, &entity);
696: GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6);
697: GmshReadSize(gmsh, &numTags, 1);
698: GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
699: GmshReadInt(gmsh, tags, numTags);
700: entity->numTags = PetscMin(numTags, 4);
701: for (i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i];
702: if (dim == 0) continue;
703: GmshReadSize(gmsh, &numTags, 1);
704: GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
705: GmshReadInt(gmsh, tags, numTags);
706: }
707: }
708: return(0);
709: }
711: /*
712: $Nodes
713: numEntityBlocks(size_t) numNodes(size_t)
714: minNodeTag(size_t) maxNodeTag(size_t)
715: entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesBlock(size_t)
716: nodeTag(size_t)
717: ...
718: x(double) y(double) z(double)
719: < u(double; if parametric and entityDim = 1 or entityDim = 2) >
720: < v(double; if parametric and entityDim = 2) >
721: ...
722: ...
723: $EndNodes
724: */
725: static PetscErrorCode DMPlexCreateGmsh_ReadNodes_v41(GmshFile *gmsh, int shift, PetscInt *numVertices, double **gmsh_nodes)
726: {
727: int info[3];
728: PetscInt sizes[4], numEntityBlocks, numNodes, numNodesBlock = 0, *nodeTag = NULL, block, node, i;
729: double *coordinates;
733: GmshReadSize(gmsh, sizes, 4);
734: numEntityBlocks = sizes[0]; numNodes = sizes[1];
735: PetscMalloc1(numNodes*3, &coordinates);
736: *numVertices = numNodes;
737: *gmsh_nodes = coordinates;
738: for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) {
739: GmshReadInt(gmsh, info, 3);
740: GmshReadSize(gmsh, &numNodesBlock, 1);
741: GmshBufferGet(gmsh, numNodesBlock, sizeof(PetscInt), &nodeTag);
742: GmshReadSize(gmsh, nodeTag, numNodesBlock);
743: 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);
744: if (info[2] != 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parametric coordinates not supported");
745: GmshReadDouble(gmsh, coordinates+node*3, numNodesBlock*3);
746: }
747: return(0);
748: }
750: /*
751: $Elements
752: numEntityBlocks(size_t) numElements(size_t)
753: minElementTag(size_t) maxElementTag(size_t)
754: entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t)
755: elementTag(size_t) nodeTag(size_t) ...
756: ...
757: ...
758: $EndElements
759: */
760: static PetscErrorCode DMPlexCreateGmsh_ReadElements_v41(GmshFile *gmsh, PETSC_UNUSED int shift, GmshEntities *entities, PetscInt *numCells, GmshElement **gmsh_elems)
761: {
762: int info[3], eid, dim, cellType, *tags;
763: PetscInt sizes[4], *ibuf = NULL, numEntityBlocks, numElements, numBlockElements, numNodes, numTags, block, elem, c, p;
764: GmshEntity *entity = NULL;
765: GmshElement *elements;
769: GmshReadSize(gmsh, sizes, 4);
770: numEntityBlocks = sizes[0]; numElements = sizes[1];
771: PetscCalloc1(numElements, &elements);
772: *numCells = numElements;
773: *gmsh_elems = elements;
774: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
775: GmshReadInt(gmsh, info, 3);
776: dim = info[0]; eid = info[1]; cellType = info[2];
777: GmshEntitiesGet(entities, dim, eid, &entity);
778: numTags = entity->numTags;
779: tags = entity->tags;
780: switch (cellType) {
781: case 1: /* 2-node line */
782: numNodes = 2;
783: break;
784: case 2: /* 3-node triangle */
785: numNodes = 3;
786: break;
787: case 3: /* 4-node quadrangle */
788: numNodes = 4;
789: break;
790: case 4: /* 4-node tetrahedron */
791: numNodes = 4;
792: break;
793: case 5: /* 8-node hexahedron */
794: numNodes = 8;
795: break;
796: case 6: /* 6-node wedge */
797: numNodes = 6;
798: break;
799: case 15: /* 1-node vertex */
800: numNodes = 1;
801: break;
802: default:
803: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
804: }
805: GmshReadSize(gmsh, &numBlockElements, 1);
806: GmshBufferGet(gmsh, (1+numNodes)*numBlockElements, sizeof(PetscInt), &ibuf);
807: GmshReadSize(gmsh, ibuf, (1+numNodes)*numBlockElements);
808: for (elem = 0; elem < numBlockElements; ++elem, ++c) {
809: GmshElement *element = elements + c;
810: PetscInt *id = ibuf + elem*(1+numNodes), *nodes = id + 1;
811: element->id = id[0];
812: element->dim = dim;
813: element->cellType = cellType;
814: element->numNodes = numNodes;
815: element->numTags = numTags;
816: for (p = 0; p < numNodes; p++) element->nodes[p] = nodes[p];
817: for (p = 0; p < numTags; p++) element->tags[p] = tags[p];
818: }
819: }
820: return(0);
821: }
823: /*
824: $Periodic
825: numPeriodicLinks(size_t)
826: entityDim(int) entityTag(int) entityTagMaster(int)
827: numAffine(size_t) value(double) ...
828: numCorrespondingNodes(size_t)
829: nodeTag(size_t) nodeTagMaster(size_t)
830: ...
831: ...
832: $EndPeriodic
833: */
834: static PetscErrorCode DMPlexCreateGmsh_ReadPeriodic_v41(GmshFile *gmsh, int shift, PetscInt slaveMap[], PetscBT bt)
835: {
836: int info[3];
837: PetscInt numPeriodicLinks, numAffine, numCorrespondingNodes, *nodeTags = NULL, link, node;
838: double dbuf[16];
842: GmshReadSize(gmsh, &numPeriodicLinks, 1);
843: for (link = 0; link < numPeriodicLinks; ++link) {
844: GmshReadInt(gmsh, info, 3);
845: GmshReadSize(gmsh, &numAffine, 1);
846: GmshReadDouble(gmsh, dbuf, numAffine);
847: GmshReadSize(gmsh, &numCorrespondingNodes, 1);
848: GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags);
849: GmshReadSize(gmsh, nodeTags, numCorrespondingNodes*2);
850: for (node = 0; node < numCorrespondingNodes; ++node) {
851: PetscInt slaveNode = nodeTags[node*2+0] - shift;
852: PetscInt masterNode = nodeTags[node*2+1] - shift;
853: slaveMap[slaveNode] = masterNode;
854: PetscBTSet(bt, slaveNode);
855: PetscBTSet(bt, masterNode);
856: }
857: }
858: return(0);
859: }
861: /*
862: $MeshFormat // same as MSH version 2
863: version(ASCII double; currently 4.1)
864: file-type(ASCII int; 0 for ASCII mode, 1 for binary mode)
865: data-size(ASCII int; sizeof(size_t))
866: < int with value one; only in binary mode, to detect endianness >
867: $EndMeshFormat
868: */
869: static PetscErrorCode DMPlexCreateGmsh_ReadMeshFormat(GmshFile *gmsh)
870: {
871: char line[PETSC_MAX_PATH_LEN];
872: int snum, fileType, fileFormat, dataSize, checkEndian;
873: float version;
877: GmshReadString(gmsh, line, 3);
878: snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
879: if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
880: 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);
881: if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f not supported", (double)version);
882: 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);
883: if (gmsh->binary && !fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Viewer is binary but Gmsh file is ASCII");
884: if (!gmsh->binary && fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Viewer is ASCII but Gmsh file is binary");
885: fileFormat = (int)(version*10); /* XXX Should use (int)roundf(version*10) ? */
886: 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);
887: 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);
888: gmsh->fileFormat = fileFormat;
889: gmsh->dataSize = dataSize;
890: gmsh->byteSwap = PETSC_FALSE;
891: if (gmsh->binary) {
892: GmshReadInt(gmsh, &checkEndian, 1);
893: if (checkEndian != 1) {
894: PetscByteSwap(&checkEndian, PETSC_ENUM, 1);
895: if (checkEndian != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to detect endianness in Gmsh file header: %s", line);
896: gmsh->byteSwap = PETSC_TRUE;
897: }
898: }
899: return(0);
900: }
902: /*
903: PhysicalNames
904: numPhysicalNames(ASCII int)
905: dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max)
906: ...
907: $EndPhysicalNames
908: */
909: static PetscErrorCode DMPlexCreateGmsh_ReadPhysicalNames(GmshFile *gmsh)
910: {
911: char line[PETSC_MAX_PATH_LEN], name[128+2], *p, *q;
912: int snum, numRegions, region, dim, tag;
916: GmshReadString(gmsh, line, 1);
917: snum = sscanf(line, "%d", &numRegions);
918: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
919: for (region = 0; region < numRegions; ++region) {
920: GmshReadString(gmsh, line, 2);
921: snum = sscanf(line, "%d %d", &dim, &tag);
922: if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
923: GmshReadString(gmsh, line, -(PetscInt)sizeof(line));
924: PetscStrchr(line, '"', &p);
925: if (!p) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
926: PetscStrrchr(line, '"', &q);
927: if (q == p) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
928: PetscStrncpy(name, p+1, (size_t)(q-p-1));
929: }
930: return(0);
931: }
933: /*@C
934: DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file
936: + comm - The MPI communicator
937: . filename - Name of the Gmsh file
938: - interpolate - Create faces and edges in the mesh
940: Output Parameter:
941: . dm - The DM object representing the mesh
943: Level: beginner
945: .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
946: @*/
947: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
948: {
949: PetscViewer viewer;
950: PetscMPIInt rank;
951: int fileType;
952: PetscViewerType vtype;
953: PetscErrorCode ierr;
956: MPI_Comm_rank(comm, &rank);
958: /* Determine Gmsh file type (ASCII or binary) from file header */
959: if (!rank) {
960: GmshFile gmsh_, *gmsh = &gmsh_;
961: char line[PETSC_MAX_PATH_LEN];
962: int snum;
963: float version;
965: PetscMemzero(gmsh,sizeof(GmshFile));
966: PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer);
967: PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII);
968: PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ);
969: PetscViewerFileSetName(gmsh->viewer, filename);
970: /* Read only the first two lines of the Gmsh file */
971: GmshReadSection(gmsh, line);
972: GmshExpect(gmsh, "$MeshFormat", line);
973: GmshReadString(gmsh, line, 2);
974: snum = sscanf(line, "%f %d", &version, &fileType);
975: if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
976: 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);
977: if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f not supported", (double)version);
978: 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);
979: PetscViewerDestroy(&gmsh->viewer);
980: }
981: MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);
982: vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;
984: /* Create appropriate viewer and build plex */
985: PetscViewerCreate(comm, &viewer);
986: PetscViewerSetType(viewer, vtype);
987: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
988: PetscViewerFileSetName(viewer, filename);
989: DMPlexCreateGmsh(comm, viewer, interpolate, dm);
990: PetscViewerDestroy(&viewer);
991: return(0);
992: }
994: /*@
995: DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer
997: Collective on comm
999: Input Parameters:
1000: + comm - The MPI communicator
1001: . viewer - The Viewer associated with a Gmsh file
1002: - interpolate - Create faces and edges in the mesh
1004: Output Parameter:
1005: . dm - The DM object representing the mesh
1007: Note: http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
1009: Level: beginner
1011: .keywords: mesh,Gmsh
1012: .seealso: DMPLEX, DMCreate()
1013: @*/
1014: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1015: {
1016: PetscViewer parentviewer = NULL;
1017: double *coordsIn = NULL;
1018: GmshEntities *entities = NULL;
1019: GmshElement *gmsh_elem = NULL;
1020: PetscSection coordSection;
1021: Vec coordinates;
1022: PetscBT periodicV = NULL, periodicC = NULL;
1023: PetscScalar *coords;
1024: PetscInt dim = 0, embedDim = -1, coordSize, c, v, d, cell, *periodicMap = NULL, *periodicMapI = NULL, *hybridMap = NULL, cMax = PETSC_DETERMINE;
1025: PetscInt numVertices = 0, numCells = 0, trueNumCells = 0;
1026: int i, shift = 1;
1027: PetscMPIInt rank;
1028: PetscBool binary, zerobase = PETSC_FALSE, periodic = PETSC_FALSE, usemarker = PETSC_FALSE;
1029: PetscBool enable_hybrid = PETSC_FALSE;
1033: MPI_Comm_rank(comm, &rank);
1034: PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_hybrid", &enable_hybrid, NULL);
1035: PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_periodic", &periodic, NULL);
1036: PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_use_marker", &usemarker, NULL);
1037: PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_zero_base", &zerobase, NULL);
1038: PetscOptionsGetInt(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_spacedim", &embedDim, NULL);
1039: if (zerobase) shift = 0;
1041: DMCreate(comm, dm);
1042: DMSetType(*dm, DMPLEX);
1043: PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);
1045: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);
1047: /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1048: if (binary) {
1049: parentviewer = viewer;
1050: PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1051: }
1053: if (!rank) {
1054: GmshFile gmsh_, *gmsh = &gmsh_;
1055: char line[PETSC_MAX_PATH_LEN];
1056: PetscBool match;
1058: PetscMemzero(gmsh,sizeof(GmshFile));
1059: gmsh->viewer = viewer;
1060: gmsh->binary = binary;
1062: /* Read mesh format */
1063: GmshReadSection(gmsh, line);
1064: GmshExpect(gmsh, "$MeshFormat", line);
1065: DMPlexCreateGmsh_ReadMeshFormat(gmsh);
1066: GmshReadEndSection(gmsh, "$EndMeshFormat", line);
1068: /* OPTIONAL Read physical names */
1069: GmshReadSection(gmsh, line);
1070: GmshMatch(gmsh,"$PhysicalNames", line, &match);
1071: if (match) {
1072: DMPlexCreateGmsh_ReadPhysicalNames(gmsh);
1073: GmshReadEndSection(gmsh, "$EndPhysicalNames", line);
1074: /* Initial read for entity section */
1075: GmshReadSection(gmsh, line);
1076: }
1078: /* Read entities */
1079: if (gmsh->fileFormat >= 40) {
1080: GmshExpect(gmsh, "$Entities", line);
1081: switch (gmsh->fileFormat) {
1082: case 41: DMPlexCreateGmsh_ReadEntities_v41(gmsh, &entities); break;
1083: default: DMPlexCreateGmsh_ReadEntities_v40(gmsh, &entities); break;
1084: }
1085: GmshReadEndSection(gmsh, "$EndEntities", line);
1086: /* Initial read for nodes section */
1087: GmshReadSection(gmsh, line);
1088: }
1090: /* Read nodes */
1091: GmshExpect(gmsh, "$Nodes", line);
1092: switch (gmsh->fileFormat) {
1093: case 41: DMPlexCreateGmsh_ReadNodes_v41(gmsh, shift, &numVertices, &coordsIn); break;
1094: case 40: DMPlexCreateGmsh_ReadNodes_v40(gmsh, shift, &numVertices, &coordsIn); break;
1095: default: DMPlexCreateGmsh_ReadNodes_v22(gmsh, shift, &numVertices, &coordsIn); break;
1096: }
1097: GmshReadEndSection(gmsh, "$EndNodes", line);
1099: /* Read elements */
1100: GmshReadSection(gmsh, line);;
1101: GmshExpect(gmsh, "$Elements", line);
1102: switch (gmsh->fileFormat) {
1103: case 41: DMPlexCreateGmsh_ReadElements_v41(gmsh, shift, entities, &numCells, &gmsh_elem); break;
1104: case 40: DMPlexCreateGmsh_ReadElements_v40(gmsh, shift, entities, &numCells, &gmsh_elem); break;
1105: default: DMPlexCreateGmsh_ReadElements_v22(gmsh, shift, &numCells, &gmsh_elem); break;
1106: }
1107: GmshReadEndSection(gmsh, "$EndElements", line);
1109: /* OPTIONAL Read periodic section */
1110: if (periodic) {
1111: PetscInt pVert, *periodicMapT, *aux;
1113: PetscMalloc1(numVertices, &periodicMapT);
1114: PetscBTCreate(numVertices, &periodicV);
1115: for (i = 0; i < numVertices; i++) periodicMapT[i] = i;
1117: GmshReadSection(gmsh, line);
1118: GmshExpect(gmsh, "$Periodic", line);
1119: switch (gmsh->fileFormat) {
1120: case 41: DMPlexCreateGmsh_ReadPeriodic_v41(gmsh, shift, periodicMapT, periodicV); break;
1121: default: DMPlexCreateGmsh_ReadPeriodic_v40(gmsh, shift, periodicMapT, periodicV); break;
1122: }
1123: GmshReadEndSection(gmsh, "$EndPeriodic", line);
1125: /* we may have slaves of slaves */
1126: for (i = 0; i < numVertices; i++) {
1127: while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) {
1128: periodicMapT[i] = periodicMapT[periodicMapT[i]];
1129: }
1130: }
1131: /* periodicMap : from old to new numbering (periodic vertices excluded)
1132: periodicMapI: from new to old numbering */
1133: PetscMalloc1(numVertices, &periodicMap);
1134: PetscMalloc1(numVertices, &periodicMapI);
1135: PetscMalloc1(numVertices, &aux);
1136: for (i = 0, pVert = 0; i < numVertices; i++) {
1137: if (periodicMapT[i] != i) {
1138: pVert++;
1139: } else {
1140: aux[i] = i - pVert;
1141: periodicMapI[i - pVert] = i;
1142: }
1143: }
1144: for (i = 0 ; i < numVertices; i++) {
1145: periodicMap[i] = aux[periodicMapT[i]];
1146: }
1147: PetscFree(periodicMapT);
1148: PetscFree(aux);
1149: /* remove periodic vertices */
1150: numVertices = numVertices - pVert;
1151: }
1153: GmshEntitiesDestroy(&entities);
1154: PetscFree(gmsh->wbuf);
1155: PetscFree(gmsh->sbuf);
1156: }
1158: if (parentviewer) {
1159: PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1160: }
1162: if (!rank) {
1163: PetscBool hybrid = PETSC_FALSE;
1165: for (trueNumCells = 0, c = 0; c < numCells; ++c) {
1166: int on = -1;
1167: if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;}
1168: if (gmsh_elem[c].dim == dim) {hybrid = (trueNumCells ? (on != gmsh_elem[c].cellType ? on = gmsh_elem[c].cellType,PETSC_TRUE : hybrid) : (on = gmsh_elem[c].cellType, PETSC_FALSE) ); trueNumCells++;}
1169: /* wedges always indicate an hybrid mesh in PLEX */
1170: if (on == 6 || on == 13) hybrid = PETSC_TRUE;
1171: }
1172: /* Renumber cells for hybrid grids */
1173: if (hybrid && enable_hybrid) {
1174: PetscInt hc1 = 0, hc2 = 0, *hybridCells1 = NULL, *hybridCells2 = NULL;
1175: PetscInt cell, tn, *tp;
1176: int n1 = 0,n2 = 0;
1178: PetscMalloc1(trueNumCells, &hybridCells1);
1179: PetscMalloc1(trueNumCells, &hybridCells2);
1180: for (cell = 0, c = 0; c < numCells; ++c) {
1181: if (gmsh_elem[c].dim == dim) {
1182: if (!n1) n1 = gmsh_elem[c].cellType;
1183: else if (!n2 && n1 != gmsh_elem[c].cellType) n2 = gmsh_elem[c].cellType;
1185: if (gmsh_elem[c].cellType == n1) hybridCells1[hc1++] = cell;
1186: else if (gmsh_elem[c].cellType == n2) hybridCells2[hc2++] = cell;
1187: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle more than 2 cell types");
1188: cell++;
1189: }
1190: }
1192: switch (n1) {
1193: case 2: /* triangles */
1194: case 9:
1195: switch (n2) {
1196: case 0: /* single type mesh */
1197: case 3: /* quads */
1198: case 10:
1199: break;
1200: default:
1201: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1202: }
1203: break;
1204: case 3:
1205: case 10:
1206: switch (n2) {
1207: case 0: /* single type mesh */
1208: case 2: /* swap since we list simplices first */
1209: case 9:
1210: tn = hc1;
1211: hc1 = hc2;
1212: hc2 = tn;
1214: tp = hybridCells1;
1215: hybridCells1 = hybridCells2;
1216: hybridCells2 = tp;
1217: break;
1218: default:
1219: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1220: }
1221: break;
1222: case 4: /* tetrahedra */
1223: case 11:
1224: switch (n2) {
1225: case 0: /* single type mesh */
1226: case 6: /* wedges */
1227: case 13:
1228: break;
1229: default:
1230: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1231: }
1232: break;
1233: case 6:
1234: case 13:
1235: switch (n2) {
1236: case 0: /* single type mesh */
1237: case 4: /* swap since we list simplices first */
1238: case 11:
1239: tn = hc1;
1240: hc1 = hc2;
1241: hc2 = tn;
1243: tp = hybridCells1;
1244: hybridCells1 = hybridCells2;
1245: hybridCells2 = tp;
1246: break;
1247: default:
1248: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1249: }
1250: break;
1251: default:
1252: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1253: }
1254: cMax = hc1;
1255: PetscMalloc1(trueNumCells, &hybridMap);
1256: for (cell = 0; cell < hc1; cell++) hybridMap[hybridCells1[cell]] = cell;
1257: for (cell = 0; cell < hc2; cell++) hybridMap[hybridCells2[cell]] = cell + hc1;
1258: PetscFree(hybridCells1);
1259: PetscFree(hybridCells2);
1260: }
1262: }
1264: /* Allocate the cell-vertex mesh */
1265: DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
1266: for (cell = 0, c = 0; c < numCells; ++c) {
1267: if (gmsh_elem[c].dim == dim) {
1268: DMPlexSetConeSize(*dm, hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].numNodes);
1269: cell++;
1270: }
1271: }
1272: DMSetUp(*dm);
1273: /* Add cell-vertex connections */
1274: for (cell = 0, c = 0; c < numCells; ++c) {
1275: if (gmsh_elem[c].dim == dim) {
1276: PetscInt pcone[8], corner;
1277: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1278: const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
1279: pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells;
1280: }
1281: if (dim == 3) {
1282: /* Tetrahedra are inverted */
1283: if (gmsh_elem[c].cellType == 4) {
1284: PetscInt tmp = pcone[0];
1285: pcone[0] = pcone[1];
1286: pcone[1] = tmp;
1287: }
1288: /* Hexahedra are inverted */
1289: if (gmsh_elem[c].cellType == 5) {
1290: PetscInt tmp = pcone[1];
1291: pcone[1] = pcone[3];
1292: pcone[3] = tmp;
1293: }
1294: /* Prisms are inverted */
1295: if (gmsh_elem[c].cellType == 6) {
1296: PetscInt tmp;
1298: tmp = pcone[1];
1299: pcone[1] = pcone[2];
1300: pcone[2] = tmp;
1301: tmp = pcone[4];
1302: pcone[4] = pcone[5];
1303: pcone[5] = tmp;
1304: }
1305: } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */
1306: /* quads are input to PLEX as prisms */
1307: if (gmsh_elem[c].cellType == 3) {
1308: PetscInt tmp = pcone[2];
1309: pcone[2] = pcone[3];
1310: pcone[3] = tmp;
1311: }
1312: }
1313: DMPlexSetCone(*dm, hybridMap ? hybridMap[cell] : cell, pcone);
1314: cell++;
1315: }
1316: }
1317: MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
1318: DMSetDimension(*dm, dim);
1319: DMPlexSetHybridBounds(*dm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
1320: DMPlexSymmetrize(*dm);
1321: DMPlexStratify(*dm);
1322: if (interpolate) {
1323: DM idm;
1325: DMPlexInterpolate(*dm, &idm);
1326: DMDestroy(dm);
1327: *dm = idm;
1328: }
1330: if (usemarker && !interpolate && dim > 1) SETERRQ(comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation");
1331: if (!rank && usemarker) {
1332: PetscInt f, fStart, fEnd;
1334: DMCreateLabel(*dm, "marker");
1335: DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
1336: for (f = fStart; f < fEnd; ++f) {
1337: PetscInt suppSize;
1339: DMPlexGetSupportSize(*dm, f, &suppSize);
1340: if (suppSize == 1) {
1341: PetscInt *cone = NULL, coneSize, p;
1343: DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1344: for (p = 0; p < coneSize; p += 2) {
1345: DMSetLabelValue(*dm, "marker", cone[p], 1);
1346: }
1347: DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1348: }
1349: }
1350: }
1352: if (!rank) {
1353: PetscInt vStart, vEnd;
1355: DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
1356: for (cell = 0, c = 0; c < numCells; ++c) {
1358: /* Create face sets */
1359: if (interpolate && gmsh_elem[c].dim == dim-1) {
1360: const PetscInt *join;
1361: PetscInt joinSize, pcone[8], corner;
1362: /* Find the relevant facet with vertex joins */
1363: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1364: const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
1365: pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart;
1366: }
1367: DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, pcone, &joinSize, &join);
1368: 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);
1369: DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
1370: DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
1371: }
1373: /* Create cell sets */
1374: if (gmsh_elem[c].dim == dim) {
1375: if (gmsh_elem[c].numTags > 0) {
1376: DMSetLabelValue(*dm, "Cell Sets", hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].tags[0]);
1377: }
1378: cell++;
1379: }
1381: /* Create vertex sets */
1382: if (gmsh_elem[c].dim == 0) {
1383: if (gmsh_elem[c].numTags > 0) {
1384: const PetscInt cc = gmsh_elem[c].nodes[0] - shift;
1385: const PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart;
1386: DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);
1387: }
1388: }
1389: }
1390: }
1392: /* Create coordinates */
1393: if (embedDim < 0) embedDim = dim;
1394: DMSetCoordinateDim(*dm, embedDim);
1395: DMGetCoordinateSection(*dm, &coordSection);
1396: PetscSectionSetNumFields(coordSection, 1);
1397: PetscSectionSetFieldComponents(coordSection, 0, embedDim);
1398: if (periodicMap) { /* we need to localize coordinates on cells */
1399: PetscSectionSetChart(coordSection, 0, trueNumCells + numVertices);
1400: } else {
1401: PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);
1402: }
1403: for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
1404: PetscSectionSetDof(coordSection, v, embedDim);
1405: PetscSectionSetFieldDof(coordSection, v, 0, embedDim);
1406: }
1407: if (periodicMap) {
1408: PetscBTCreate(trueNumCells, &periodicC);
1409: for (cell = 0, c = 0; c < numCells; ++c) {
1410: if (gmsh_elem[c].dim == dim) {
1411: PetscInt corner;
1412: PetscBool pc = PETSC_FALSE;
1413: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1414: pc = (PetscBool)(pc || PetscBTLookup(periodicV, gmsh_elem[c].nodes[corner] - shift));
1415: }
1416: if (pc) {
1417: PetscInt dof = gmsh_elem[c].numNodes*embedDim;
1418: PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1419: PetscBTSet(periodicC, ucell);
1420: PetscSectionSetDof(coordSection, ucell, dof);
1421: PetscSectionSetFieldDof(coordSection, ucell, 0, dof);
1422: }
1423: cell++;
1424: }
1425: }
1426: }
1427: PetscSectionSetUp(coordSection);
1428: PetscSectionGetStorageSize(coordSection, &coordSize);
1429: VecCreate(PETSC_COMM_SELF, &coordinates);
1430: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1431: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1432: VecSetBlockSize(coordinates, embedDim);
1433: VecSetType(coordinates, VECSTANDARD);
1434: VecGetArray(coordinates, &coords);
1435: if (periodicMap) {
1436: PetscInt off;
1438: for (cell = 0, c = 0; c < numCells; ++c) {
1439: PetscInt pcone[8], corner;
1440: if (gmsh_elem[c].dim == dim) {
1441: PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1442: if (PetscUnlikely(PetscBTLookup(periodicC, ucell))) {
1443: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1444: pcone[corner] = gmsh_elem[c].nodes[corner] - shift;
1445: }
1446: if (dim == 3) {
1447: /* Tetrahedra are inverted */
1448: if (gmsh_elem[c].cellType == 4) {
1449: PetscInt tmp = pcone[0];
1450: pcone[0] = pcone[1];
1451: pcone[1] = tmp;
1452: }
1453: /* Hexahedra are inverted */
1454: if (gmsh_elem[c].cellType == 5) {
1455: PetscInt tmp = pcone[1];
1456: pcone[1] = pcone[3];
1457: pcone[3] = tmp;
1458: }
1459: /* Prisms are inverted */
1460: if (gmsh_elem[c].cellType == 6) {
1461: PetscInt tmp;
1463: tmp = pcone[1];
1464: pcone[1] = pcone[2];
1465: pcone[2] = tmp;
1466: tmp = pcone[4];
1467: pcone[4] = pcone[5];
1468: pcone[5] = tmp;
1469: }
1470: } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */
1471: /* quads are input to PLEX as prisms */
1472: if (gmsh_elem[c].cellType == 3) {
1473: PetscInt tmp = pcone[2];
1474: pcone[2] = pcone[3];
1475: pcone[3] = tmp;
1476: }
1477: }
1478: PetscSectionGetOffset(coordSection, ucell, &off);
1479: for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1480: v = pcone[corner];
1481: for (d = 0; d < embedDim; ++d) {
1482: coords[off++] = (PetscReal) coordsIn[v*3+d];
1483: }
1484: }
1485: }
1486: cell++;
1487: }
1488: }
1489: for (v = 0; v < numVertices; ++v) {
1490: PetscSectionGetOffset(coordSection, v + trueNumCells, &off);
1491: for (d = 0; d < embedDim; ++d) {
1492: coords[off+d] = (PetscReal) coordsIn[periodicMapI[v]*3+d];
1493: }
1494: }
1495: } else {
1496: for (v = 0; v < numVertices; ++v) {
1497: for (d = 0; d < embedDim; ++d) {
1498: coords[v*embedDim+d] = (PetscReal) coordsIn[v*3+d];
1499: }
1500: }
1501: }
1502: VecRestoreArray(coordinates, &coords);
1503: DMSetCoordinatesLocal(*dm, coordinates);
1504: DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);
1506: PetscFree(coordsIn);
1507: PetscFree(gmsh_elem);
1508: PetscFree(hybridMap);
1509: PetscFree(periodicMap);
1510: PetscFree(periodicMapI);
1511: PetscBTDestroy(&periodicV);
1512: PetscBTDestroy(&periodicC);
1513: VecDestroy(&coordinates);
1515: PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
1516: return(0);
1517: }