Actual source code: plexgmsh.c

petsc-master 2020-04-07
Report Typos and Errors
  1: #define PETSCDM_DLL
  2:  #include <petsc/private/dmpleximpl.h>
  3:  #include <petsc/private/hashmapi.h>

  5: typedef struct {
  6:   PetscViewer    viewer;
  7:   int            fileFormat;
  8:   int            dataSize;
  9:   PetscBool      binary;
 10:   PetscBool      byteSwap;
 11:   size_t         wlen;
 12:   void           *wbuf;
 13:   size_t         slen;
 14:   void           *sbuf;
 15: } GmshFile;

 17: static PetscErrorCode GmshBufferGet(GmshFile *gmsh, size_t count, size_t eltsize, void *buf)
 18: {
 19:   size_t         size = count * eltsize;

 23:   if (gmsh->wlen < size) {
 24:     PetscFree(gmsh->wbuf);
 25:     PetscMalloc(size, &gmsh->wbuf);
 26:     gmsh->wlen = size;
 27:   }
 28:   *(void**)buf = size ? gmsh->wbuf : NULL;
 29:   return(0);
 30: }

 32: static PetscErrorCode GmshBufferSizeGet(GmshFile *gmsh, size_t count, void *buf)
 33: {
 34:   size_t         dataSize = (size_t)gmsh->dataSize;
 35:   size_t         size = count * dataSize;

 39:   if (gmsh->slen < size) {
 40:     PetscFree(gmsh->sbuf);
 41:     PetscMalloc(size, &gmsh->sbuf);
 42:     gmsh->slen = size;
 43:   }
 44:   *(void**)buf = size ? gmsh->sbuf : NULL;
 45:   return(0);
 46: }

 48: static PetscErrorCode GmshRead(GmshFile *gmsh, void *buf, PetscInt count, PetscDataType dtype)
 49: {
 52:   PetscViewerRead(gmsh->viewer, buf, count, NULL, dtype);
 53:   if (gmsh->byteSwap) {PetscByteSwap(buf, dtype, count);}
 54:   return(0);
 55: }

 57: static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscInt count)
 58: {
 61:   PetscViewerRead(gmsh->viewer, buf, count, NULL, PETSC_STRING);
 62:   return(0);
 63: }

 65: static PetscErrorCode GmshMatch(PETSC_UNUSED GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN], PetscBool *match)
 66: {
 69:   PetscStrcmp(line, Section, match);
 70:   return(0);
 71: }

 73: static PetscErrorCode GmshExpect(GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN])
 74: {
 75:   PetscBool      match;

 79:   GmshMatch(gmsh, Section, line, &match);
 80:   if (!match) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file, expecting %s",Section);
 81:   return(0);
 82: }

 84: static PetscErrorCode GmshReadSection(GmshFile *gmsh, char line[PETSC_MAX_PATH_LEN])
 85: {
 86:   PetscBool      match;

 90:   while (PETSC_TRUE) {
 91:     GmshReadString(gmsh, line, 1);
 92:     GmshMatch(gmsh, "$Comments", line, &match);
 93:     if (!match) break;
 94:     while (PETSC_TRUE) {
 95:       GmshReadString(gmsh, line, 1);
 96:       GmshMatch(gmsh, "$EndComments", line, &match);
 97:       if (match) break;
 98:     }
 99:   }
100:   return(0);
101: }

103: static PetscErrorCode GmshReadEndSection(GmshFile *gmsh, const char EndSection[], char line[PETSC_MAX_PATH_LEN])
104: {
107:   GmshReadString(gmsh, line, 1);
108:   GmshExpect(gmsh, EndSection, line);
109:   return(0);
110: }

112: static PetscErrorCode GmshReadSize(GmshFile *gmsh, PetscInt *buf, PetscInt count)
113: {
114:   PetscInt       i;
115:   size_t         dataSize = (size_t)gmsh->dataSize;

119:   if (dataSize == sizeof(PetscInt)) {
120:     GmshRead(gmsh, buf, count, PETSC_INT);
121:   } else  if (dataSize == sizeof(int)) {
122:     int *ibuf = NULL;
123:     GmshBufferSizeGet(gmsh, count, &ibuf);
124:     GmshRead(gmsh, ibuf, count, PETSC_ENUM);
125:     for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
126:   } else  if (dataSize == sizeof(long)) {
127:     long *ibuf = NULL;
128:     GmshBufferSizeGet(gmsh, count, &ibuf);
129:     GmshRead(gmsh, ibuf, count, PETSC_LONG);
130:     for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
131:   } else if (dataSize == sizeof(PetscInt64)) {
132:     PetscInt64 *ibuf = NULL;
133:     GmshBufferSizeGet(gmsh, count, &ibuf);
134:     GmshRead(gmsh, ibuf, count, PETSC_INT64);
135:     for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
136:   }
137:   return(0);
138: }

140: static PetscErrorCode GmshReadInt(GmshFile *gmsh, int *buf, PetscInt count)
141: {
144:   GmshRead(gmsh, buf, count, PETSC_ENUM);
145:   return(0);
146: }

148: static PetscErrorCode GmshReadDouble(GmshFile *gmsh, double *buf, PetscInt count)
149: {
152:   GmshRead(gmsh, buf, count, PETSC_DOUBLE);
153:   return(0);
154: }

156: typedef struct {
157:   PetscInt id;       /* Entity number */
158:   PetscInt dim;      /* Entity dimension */
159:   double   bbox[6];  /* Bounding box */
160:   PetscInt numTags;  /* Size of tag array */
161:   int      tags[4];  /* Tag array */
162: } GmshEntity;

164: typedef struct {
165:   GmshEntity *entity[4];
166:   PetscHMapI  entityMap[4];
167: } GmshEntities;

169: static PetscErrorCode GmshEntitiesCreate(PetscInt count[4], GmshEntities **entities)
170: {
171:   PetscInt       dim;

175:   PetscNew(entities);
176:   for (dim = 0; dim < 4; ++dim) {
177:     PetscCalloc1(count[dim], &(*entities)->entity[dim]);
178:     PetscHMapICreate(&(*entities)->entityMap[dim]);
179:   }
180:   return(0);
181: }

183: static PetscErrorCode GmshEntitiesAdd(GmshEntities *entities, PetscInt index, PetscInt dim, PetscInt eid, GmshEntity** entity)
184: {
187:   PetscHMapISet(entities->entityMap[dim], eid, index);
188:   entities->entity[dim][index].dim = dim;
189:   entities->entity[dim][index].id  = eid;
190:   if (entity) *entity = &entities->entity[dim][index];
191:   return(0);
192: }

194: static PetscErrorCode GmshEntitiesGet(GmshEntities *entities, PetscInt dim, PetscInt eid, GmshEntity** entity)
195: {
196:   PetscInt       index;

200:   PetscHMapIGet(entities->entityMap[dim], eid, &index);
201:   *entity = &entities->entity[dim][index];
202:   return(0);
203: }

205: static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities)
206: {
207:   PetscInt       dim;

211:   if (!*entities) return(0);
212:   for (dim = 0; dim < 4; ++dim) {
213:     PetscFree((*entities)->entity[dim]);
214:     PetscHMapIDestroy(&(*entities)->entityMap[dim]);
215:   }
216:   PetscFree((*entities));
217:   return(0);
218: }

220: typedef struct {
221:   PetscInt id;       /* Entity number */
222:   PetscInt dim;      /* Entity dimension */
223:   PetscInt cellType; /* Cell type */
224:   PetscInt numNodes; /* Size of node array */
225:   PetscInt nodes[8]; /* Node array */
226:   PetscInt numTags;  /* Size of tag array */
227:   int      tags[4];  /* Tag array */
228: } GmshElement;

230: static PetscErrorCode DMPlexCreateGmsh_ReadNodes_v22(GmshFile *gmsh, int shift, PetscInt *numVertices, double **gmsh_nodes)
231: {
232:   PetscViewer    viewer = gmsh->viewer;
233:   PetscBool      byteSwap = gmsh->byteSwap;
234:   char           line[PETSC_MAX_PATH_LEN];
235:   int            v, num, nid, snum;
236:   double         *coordinates;

240:   PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
241:   snum = sscanf(line, "%d", &num);
242:   if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
243:   PetscMalloc1(num*3, &coordinates);
244:   *numVertices = num;
245:   *gmsh_nodes = coordinates;
246:   for (v = 0; v < num; ++v) {
247:     double *xyz = coordinates + v*3;
248:     PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
249:     if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
250:     if (nid != v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %d", nid, v+shift);
251:     PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
252:     if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
253:   }
254:   return(0);
255: }

257: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
258:    file contents multiple times to figure out the true number of cells and facets
259:    in the given mesh. To make this more efficient we read the file contents only
260:    once and store them in memory, while determining the true number of cells. */
261: static PetscErrorCode DMPlexCreateGmsh_ReadElements_v22(GmshFile* gmsh, PETSC_UNUSED int shift, PetscInt *numCells, GmshElement **gmsh_elems)
262: {
263:   PetscViewer    viewer = gmsh->viewer;
264:   PetscBool      binary = gmsh->binary;
265:   PetscBool      byteSwap = gmsh->byteSwap;
266:   char           line[PETSC_MAX_PATH_LEN];
267:   GmshElement   *elements;
268:   int            i, c, p, num, ibuf[1+4+512], snum;
269:   int            cellType, dim, numNodes, numNodesIgnore, numElem, numTags;

273:   PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
274:   snum = sscanf(line, "%d", &num);
275:   if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
276:   PetscMalloc1(num, &elements);
277:   *numCells = num;
278:   *gmsh_elems = elements;
279:   for (c = 0; c < num;) {
280:     PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
281:     if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 3);}
282:     if (binary) {
283:       cellType = ibuf[0];
284:       numElem = ibuf[1];
285:       numTags = ibuf[2];
286:     } else {
287:       elements[c].id = ibuf[0];
288:       cellType = ibuf[1];
289:       numTags = ibuf[2];
290:       numElem = 1;
291:     }
292:     /* http://gmsh.info/doc/texinfo/gmsh.html#MSH-ASCII-file-format */
293:     numNodesIgnore = 0;
294:     switch (cellType) {
295:     case 1: /* 2-node line */
296:       dim = 1;
297:       numNodes = 2;
298:       break;
299:     case 2: /* 3-node triangle */
300:       dim = 2;
301:       numNodes = 3;
302:       break;
303:     case 3: /* 4-node quadrangle */
304:       dim = 2;
305:       numNodes = 4;
306:       break;
307:     case 4: /* 4-node tetrahedron */
308:       dim  = 3;
309:       numNodes = 4;
310:       break;
311:     case 5: /* 8-node hexahedron */
312:       dim = 3;
313:       numNodes = 8;
314:       break;
315:     case 6: /* 6-node wedge */
316:       dim = 3;
317:       numNodes = 6;
318:       break;
319:     case 8: /* 3-node 2nd order line */
320:       dim = 1;
321:       numNodes = 2;
322:       numNodesIgnore = 1;
323:       break;
324:     case 9: /* 6-node 2nd order triangle */
325:       dim = 2;
326:       numNodes = 3;
327:       numNodesIgnore = 3;
328:       break;
329:     case 10: /* 9-node 2nd order quadrangle */
330:       dim = 2;
331:       numNodes = 4;
332:       numNodesIgnore = 5;
333:       break;
334:     case 11: /* 10-node 2nd order tetrahedron */
335:       dim  = 3;
336:       numNodes = 4;
337:       numNodesIgnore = 6;
338:       break;
339:     case 12: /* 27-node 2nd order hexhedron */
340:       dim  = 3;
341:       numNodes = 8;
342:       numNodesIgnore = 19;
343:       break;
344:     case 13: /* 18-node 2nd wedge */
345:       dim = 3;
346:       numNodes = 6;
347:       numNodesIgnore = 12;
348:       break;
349:     case 15: /* 1-node vertex */
350:       dim = 0;
351:       numNodes = 1;
352:       break;
353:     case 7: /* 5-node pyramid */
354:     case 14: /* 14-node 2nd order pyramid */
355:     default:
356:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
357:     }
358:     if (binary) {
359:       const int nint = 1 + numTags + numNodes + numNodesIgnore;
360:       /* Loop over element blocks */
361:       for (i = 0; i < numElem; ++i, ++c) {
362:         PetscViewerRead(viewer, ibuf, nint, NULL, PETSC_ENUM);
363:         if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, nint);}
364:         elements[c].dim = dim;
365:         elements[c].numNodes = numNodes;
366:         elements[c].numTags = numTags;
367:         elements[c].id = ibuf[0];
368:         elements[c].cellType = cellType;
369:         for (p = 0; p < numTags;  p++) elements[c].tags[p]  = ibuf[1 + p];
370:         for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
371:       }
372:     } else {
373:       const int nint = numTags + numNodes + numNodesIgnore;
374:       elements[c].dim = dim;
375:       elements[c].numNodes = numNodes;
376:       elements[c].numTags = numTags;
377:       elements[c].cellType = cellType;
378:       PetscViewerRead(viewer, ibuf, nint, NULL, PETSC_ENUM);
379:       for (p = 0; p < numTags;  p++) elements[c].tags[p]  = ibuf[p];
380:       for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[numTags + p];
381:       c++;
382:     }
383:   }
384:   return(0);
385: }

387: /*
388: $Entities
389:   numPoints(unsigned long) numCurves(unsigned long)
390:     numSurfaces(unsigned long) numVolumes(unsigned long)
391:   // points
392:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
393:     boxMaxX(double) boxMaxY(double) boxMaxZ(double)
394:     numPhysicals(unsigned long) phyisicalTag[...](int)
395:   ...
396:   // curves
397:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
398:      boxMaxX(double) boxMaxY(double) boxMaxZ(double)
399:      numPhysicals(unsigned long) physicalTag[...](int)
400:      numBREPVert(unsigned long) tagBREPVert[...](int)
401:   ...
402:   // surfaces
403:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
404:     boxMaxX(double) boxMaxY(double) boxMaxZ(double)
405:     numPhysicals(unsigned long) physicalTag[...](int)
406:     numBREPCurve(unsigned long) tagBREPCurve[...](int)
407:   ...
408:   // volumes
409:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
410:     boxMaxX(double) boxMaxY(double) boxMaxZ(double)
411:     numPhysicals(unsigned long) physicalTag[...](int)
412:     numBREPSurfaces(unsigned long) tagBREPSurfaces[...](int)
413:   ...
414: $EndEntities
415: */
416: static PetscErrorCode DMPlexCreateGmsh_ReadEntities_v40(GmshFile *gmsh, GmshEntities **entities)
417: {
418:   PetscViewer    viewer = gmsh->viewer;
419:   PetscBool      byteSwap = gmsh->byteSwap;
420:   long           index, num, lbuf[4];
421:   int            dim, eid, numTags, *ibuf, t;
422:   PetscInt       count[4], i;
423:   GmshEntity     *entity = NULL;

427:   PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG);
428:   if (byteSwap) {PetscByteSwap(lbuf, PETSC_LONG, 4);}
429:   for (i = 0; i < 4; ++i) count[i] = lbuf[i];
430:   GmshEntitiesCreate(count, entities);
431:   for (dim = 0; dim < 4; ++dim) {
432:     for (index = 0; index < count[dim]; ++index) {
433:       PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM);
434:       if (byteSwap) {PetscByteSwap(&eid, PETSC_ENUM, 1);}
435:       GmshEntitiesAdd(*entities, (PetscInt)index, dim, eid, &entity);
436:       PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE);
437:       if (byteSwap) {PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6);}
438:       PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
439:       if (byteSwap) {PetscByteSwap(&num, PETSC_LONG, 1);}
440:       GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
441:       PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
442:       if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, num);}
443:       entity->numTags = numTags = (int) PetscMin(num, 4);
444:       for (t = 0; t < numTags; ++t) entity->tags[t] = ibuf[t];
445:       if (dim == 0) continue;
446:       PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
447:       if (byteSwap) {PetscByteSwap(&num, PETSC_LONG, 1);}
448:       GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
449:       PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
450:       if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, num);}
451:     }
452:   }
453:   return(0);
454: }

456: /*
457: $Nodes
458:   numEntityBlocks(unsigned long) numNodes(unsigned long)
459:   tagEntity(int) dimEntity(int) typeNode(int) numNodes(unsigned long)
460:     tag(int) x(double) y(double) z(double)
461:     ...
462:   ...
463: $EndNodes
464: */
465: static PetscErrorCode DMPlexCreateGmsh_ReadNodes_v40(GmshFile *gmsh, int shift, PetscInt *numVertices, double **gmsh_nodes)
466: {
467:   PetscViewer    viewer = gmsh->viewer;
468:   PetscBool      byteSwap = gmsh->byteSwap;
469:   long           block, node, v, numEntityBlocks, numTotalNodes, numNodes;
470:   int            info[3], nid;
471:   double         *coordinates;

475:   PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
476:   if (byteSwap) {PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);}
477:   PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG);
478:   if (byteSwap) {PetscByteSwap(&numTotalNodes, PETSC_LONG, 1);}
479:   PetscMalloc1(numTotalNodes*3, &coordinates);
480:   *numVertices = numTotalNodes;
481:   *gmsh_nodes = coordinates;
482:   for (v = 0, block = 0; block < numEntityBlocks; ++block) {
483:     PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
484:     PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG);
485:     if (byteSwap) {PetscByteSwap(&numNodes, PETSC_LONG, 1);}
486:     if (gmsh->binary) {
487:       size_t nbytes = sizeof(int) + 3*sizeof(double);
488:       char   *cbuf = NULL; /* dummy value to prevent warning from compiler about possible unitilized value */
489:       GmshBufferGet(gmsh, numNodes, nbytes, &cbuf);
490:       PetscViewerRead(viewer, cbuf, numNodes*nbytes, NULL, PETSC_CHAR);
491:       for (node = 0; node < numNodes; ++node, ++v) {
492:         char   *cnid = cbuf + node*nbytes, *cxyz = cnid + sizeof(int);
493:         double *xyz = coordinates + v*3;
494:         if (!PetscBinaryBigEndian()) {PetscByteSwap(cnid, PETSC_ENUM, 1);}
495:         if (!PetscBinaryBigEndian()) {PetscByteSwap(cxyz, PETSC_DOUBLE, 3);}
496:         PetscMemcpy(&nid, cnid, sizeof(int));
497:         PetscMemcpy(xyz, cxyz, 3*sizeof(double));
498:         if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
499:         if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
500:         if ((long)nid != v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %ld", nid, v+shift);
501:       }
502:     } else {
503:       for (node = 0; node < numNodes; ++node, ++v) {
504:         double *xyz = coordinates + v*3;
505:         PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
506:         if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
507:         if ((long)nid != v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %ld", nid, v+shift);
508:         PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
509:         if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
510:       }
511:     }
512:   }
513:   return(0);
514: }

516: /*
517: $Elements
518:   numEntityBlocks(unsigned long) numElements(unsigned long)
519:   tagEntity(int) dimEntity(int) typeEle(int) numElements(unsigned long)
520:     tag(int) numVert[...](int)
521:     ...
522:   ...
523: $EndElements
524: */
525: static PetscErrorCode DMPlexCreateGmsh_ReadElements_v40(GmshFile *gmsh, PETSC_UNUSED int shift, GmshEntities *entities, PetscInt *numCells, GmshElement **gmsh_elems)
526: {
527:   PetscViewer    viewer = gmsh->viewer;
528:   PetscBool      byteSwap = gmsh->byteSwap;
529:   long           c, block, numEntityBlocks, numTotalElements, elem, numElements;
530:   int            p, info[3], *ibuf = NULL;
531:   int            eid, dim, numTags, *tags, cellType, numNodes;
532:   GmshEntity     *entity = NULL;
533:   GmshElement    *elements;

537:   PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
538:   if (byteSwap) {PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);}
539:   PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG);
540:   if (byteSwap) {PetscByteSwap(&numTotalElements, PETSC_LONG, 1);}
541:   PetscCalloc1(numTotalElements, &elements);
542:   *numCells = numTotalElements;
543:   *gmsh_elems = elements;
544:   for (c = 0, block = 0; block < numEntityBlocks; ++block) {
545:     PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
546:     if (byteSwap) {PetscByteSwap(info, PETSC_ENUM, 3);}
547:     eid = info[0]; dim = info[1]; cellType = info[2];
548:     GmshEntitiesGet(entities, dim, eid, &entity);
549:     numTags = entity->numTags;
550:     tags = entity->tags;
551:     switch (cellType) {
552:     case 1: /* 2-node line */
553:       numNodes = 2;
554:       break;
555:     case 2: /* 3-node triangle */
556:       numNodes = 3;
557:       break;
558:     case 3: /* 4-node quadrangle */
559:       numNodes = 4;
560:       break;
561:     case 4: /* 4-node tetrahedron */
562:       numNodes = 4;
563:       break;
564:     case 5: /* 8-node hexahedron */
565:       numNodes = 8;
566:       break;
567:     case 6: /* 6-node wedge */
568:       numNodes = 6;
569:       break;
570:     case 15: /* 1-node vertex */
571:       numNodes = 1;
572:       break;
573:     default:
574:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
575:     }
576:     PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG);
577:     if (byteSwap) {PetscByteSwap(&numElements, PETSC_LONG, 1);}
578:     GmshBufferGet(gmsh, (1+numNodes)*numElements, sizeof(int), &ibuf);
579:     PetscViewerRead(viewer, ibuf, (1+numNodes)*numElements, NULL, PETSC_ENUM);
580:     if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, (1+numNodes)*numElements);}
581:     for (elem = 0; elem < numElements; ++elem, ++c) {
582:       int *id = ibuf + elem*(1+numNodes), *nodes = id + 1;
583:       GmshElement *element = elements + c;
584:       element->dim = dim;
585:       element->cellType = cellType;
586:       element->numNodes = numNodes;
587:       element->numTags = numTags;
588:       element->id = id[0];
589:       for (p = 0; p < numNodes; p++) element->nodes[p] = nodes[p];
590:       for (p = 0; p < numTags;  p++) element->tags[p]  = tags[p];
591:     }
592:   }
593:   return(0);
594: }

596: static PetscErrorCode DMPlexCreateGmsh_ReadPeriodic_v40(GmshFile *gmsh, int shift, PetscInt slaveMap[], PetscBT bt)
597: {
598:   PetscViewer    viewer = gmsh->viewer;
599:   int            fileFormat = gmsh->fileFormat;
600:   PetscBool      binary = gmsh->binary;
601:   PetscBool      byteSwap = gmsh->byteSwap;
602:   int            numPeriodic, snum, i;
603:   char           line[PETSC_MAX_PATH_LEN];

607:   if (fileFormat == 22 || !binary) {
608:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
609:     snum = sscanf(line, "%d", &numPeriodic);
610:     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
611:   } else {
612:     PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM);
613:     if (byteSwap) {PetscByteSwap(&numPeriodic, PETSC_ENUM, 1);}
614:   }
615:   for (i = 0; i < numPeriodic; i++) {
616:     int    ibuf[3], slaveDim = -1, slaveTag = -1, masterTag = -1, slaveNode, masterNode;
617:     long   j, nNodes;
618:     double affine[16];

620:     if (fileFormat == 22 || !binary) {
621:       PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
622:       snum = sscanf(line, "%d %d %d", &slaveDim, &slaveTag, &masterTag);
623:       if (snum != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
624:     } else {
625:       PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
626:       if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 3);}
627:       slaveDim = ibuf[0]; slaveTag = ibuf[1]; masterTag = ibuf[2];
628:     }
629:     (void)slaveDim; (void)slaveTag; (void)masterTag; /* unused */

631:     if (fileFormat == 22 || !binary) {
632:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
633:       snum = sscanf(line, "%ld", &nNodes);
634:       if (snum != 1) { /* discard transformation and try again */
635:         PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
636:         PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
637:         snum = sscanf(line, "%ld", &nNodes);
638:         if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
639:       }
640:     } else {
641:       PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
642:       if (byteSwap) {PetscByteSwap(&nNodes, PETSC_LONG, 1);}
643:       if (nNodes == -1) { /* discard transformation and try again */
644:         PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE);
645:         PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
646:         if (byteSwap) {PetscByteSwap(&nNodes, PETSC_LONG, 1);}
647:       }
648:     }

650:     for (j = 0; j < nNodes; j++) {
651:       if (fileFormat == 22 || !binary) {
652:         PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
653:         snum = sscanf(line, "%d %d", &slaveNode, &masterNode);
654:         if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
655:       } else {
656:         PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM);
657:         if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 2);}
658:         slaveNode = ibuf[0]; masterNode = ibuf[1];
659:       }
660:       slaveMap[slaveNode - shift] = masterNode - shift;
661:       PetscBTSet(bt, slaveNode - shift);
662:       PetscBTSet(bt, masterNode - shift);
663:     }
664:   }
665:   return(0);
666: }

668: /*
669: $Entities
670:   numPoints(size_t) numCurves(size_t)
671:     numSurfaces(size_t) numVolumes(size_t)
672:   pointTag(int) X(double) Y(double) Z(double)
673:     numPhysicalTags(size_t) physicalTag(int) ...
674:   ...
675:   curveTag(int) minX(double) minY(double) minZ(double)
676:     maxX(double) maxY(double) maxZ(double)
677:     numPhysicalTags(size_t) physicalTag(int) ...
678:     numBoundingPoints(size_t) pointTag(int) ...
679:   ...
680:   surfaceTag(int) minX(double) minY(double) minZ(double)
681:     maxX(double) maxY(double) maxZ(double)
682:     numPhysicalTags(size_t) physicalTag(int) ...
683:     numBoundingCurves(size_t) curveTag(int) ...
684:   ...
685:   volumeTag(int) minX(double) minY(double) minZ(double)
686:     maxX(double) maxY(double) maxZ(double)
687:     numPhysicalTags(size_t) physicalTag(int) ...
688:     numBoundngSurfaces(size_t) surfaceTag(int) ...
689:   ...
690: $EndEntities
691: */
692: static PetscErrorCode DMPlexCreateGmsh_ReadEntities_v41(GmshFile *gmsh, GmshEntities **entities)
693: {
694:   PetscInt       count[4], index, numTags, i;
695:   int            dim, eid, *tags = NULL;
696:   GmshEntity     *entity = NULL;

700:   GmshReadSize(gmsh, count, 4);
701:   GmshEntitiesCreate(count, entities);
702:   for (dim = 0; dim < 4; ++dim) {
703:     for (index = 0; index < count[dim]; ++index) {
704:       GmshReadInt(gmsh, &eid, 1);
705:       GmshEntitiesAdd(*entities, (PetscInt)index, dim, eid, &entity);
706:       GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6);
707:       GmshReadSize(gmsh, &numTags, 1);
708:       GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
709:       GmshReadInt(gmsh, tags, numTags);
710:       entity->numTags = PetscMin(numTags, 4);
711:       for (i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i];
712:       if (dim == 0) continue;
713:       GmshReadSize(gmsh, &numTags, 1);
714:       GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
715:       GmshReadInt(gmsh, tags, numTags);
716:     }
717:   }
718:   return(0);
719: }

721: /*
722: $Nodes
723:   numEntityBlocks(size_t) numNodes(size_t)
724:     minNodeTag(size_t) maxNodeTag(size_t)
725:   entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesBlock(size_t)
726:     nodeTag(size_t)
727:     ...
728:     x(double) y(double) z(double)
729:        < u(double; if parametric and entityDim = 1 or entityDim = 2) >
730:        < v(double; if parametric and entityDim = 2) >
731:     ...
732:   ...
733: $EndNodes
734: */
735: static PetscErrorCode DMPlexCreateGmsh_ReadNodes_v41(GmshFile *gmsh, int shift, PetscInt *numVertices, double **gmsh_nodes)
736: {
737:   int            info[3];
738:   PetscInt       sizes[4], numEntityBlocks, numNodes, numNodesBlock = 0, *nodeTag = NULL, block, node, i;
739:   double         *coordinates;

743:   GmshReadSize(gmsh, sizes, 4);
744:   numEntityBlocks = sizes[0]; numNodes = sizes[1];
745:   PetscMalloc1(numNodes*3, &coordinates);
746:   *numVertices = numNodes;
747:   *gmsh_nodes = coordinates;
748:   for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) {
749:     GmshReadInt(gmsh, info, 3);
750:     GmshReadSize(gmsh, &numNodesBlock, 1);
751:     GmshBufferGet(gmsh, numNodesBlock, sizeof(PetscInt), &nodeTag);
752:     GmshReadSize(gmsh, nodeTag, numNodesBlock);
753:     for (i = 0; i < numNodesBlock; ++i) if (nodeTag[i] != node+i+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %d", nodeTag[i], node+i+shift);
754:     if (info[2] != 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parametric coordinates not supported");
755:     GmshReadDouble(gmsh, coordinates+node*3, numNodesBlock*3);
756:   }
757:   return(0);
758: }

760: /*
761: $Elements
762:   numEntityBlocks(size_t) numElements(size_t)
763:     minElementTag(size_t) maxElementTag(size_t)
764:   entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t)
765:     elementTag(size_t) nodeTag(size_t) ...
766:     ...
767:   ...
768: $EndElements
769: */
770: static PetscErrorCode DMPlexCreateGmsh_ReadElements_v41(GmshFile *gmsh, PETSC_UNUSED int shift, GmshEntities *entities, PetscInt *numCells, GmshElement **gmsh_elems)
771: {
772:   int            info[3], eid, dim, cellType, *tags;
773:   PetscInt       sizes[4], *ibuf = NULL, numEntityBlocks, numElements, numBlockElements, numNodes, numTags, block, elem, c, p;
774:   GmshEntity     *entity = NULL;
775:   GmshElement    *elements;

779:   GmshReadSize(gmsh, sizes, 4);
780:   numEntityBlocks = sizes[0]; numElements = sizes[1];
781:   PetscCalloc1(numElements, &elements);
782:   *numCells = numElements;
783:   *gmsh_elems = elements;
784:   for (c = 0, block = 0; block < numEntityBlocks; ++block) {
785:     GmshReadInt(gmsh, info, 3);
786:     dim = info[0]; eid = info[1]; cellType = info[2];
787:     GmshEntitiesGet(entities, dim, eid, &entity);
788:     numTags = entity->numTags;
789:     tags = entity->tags;
790:     switch (cellType) {
791:     case 1: /* 2-node line */
792:       numNodes = 2;
793:       break;
794:     case 2: /* 3-node triangle */
795:       numNodes = 3;
796:       break;
797:     case 3: /* 4-node quadrangle */
798:       numNodes = 4;
799:       break;
800:     case 4: /* 4-node tetrahedron */
801:       numNodes = 4;
802:       break;
803:     case 5: /* 8-node hexahedron */
804:       numNodes = 8;
805:       break;
806:     case 6: /* 6-node wedge */
807:       numNodes = 6;
808:       break;
809:     case 15: /* 1-node vertex */
810:       numNodes = 1;
811:       break;
812:     default:
813:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
814:     }
815:     GmshReadSize(gmsh, &numBlockElements, 1);
816:     GmshBufferGet(gmsh, (1+numNodes)*numBlockElements, sizeof(PetscInt), &ibuf);
817:     GmshReadSize(gmsh, ibuf, (1+numNodes)*numBlockElements);
818:     for (elem = 0; elem < numBlockElements; ++elem, ++c) {
819:       GmshElement *element = elements + c;
820:       PetscInt *id = ibuf + elem*(1+numNodes), *nodes = id + 1;
821:       element->id       = id[0];
822:       element->dim      = dim;
823:       element->cellType = cellType;
824:       element->numNodes = numNodes;
825:       element->numTags  = numTags;
826:       for (p = 0; p < numNodes; p++) element->nodes[p] = nodes[p];
827:       for (p = 0; p < numTags;  p++) element->tags[p]  = tags[p];
828:     }
829:   }
830:   return(0);
831: }

833: /*
834: $Periodic
835:   numPeriodicLinks(size_t)
836:  entityDim(int) entityTag(int) entityTagMaster(int)
837:   numAffine(size_t) value(double) ...
838:   numCorrespondingNodes(size_t)
839:     nodeTag(size_t) nodeTagMaster(size_t)
840:     ...
841:   ...
842: $EndPeriodic
843: */
844: static PetscErrorCode DMPlexCreateGmsh_ReadPeriodic_v41(GmshFile *gmsh, int shift, PetscInt slaveMap[], PetscBT bt)
845: {
846:   int            info[3];
847:   PetscInt       numPeriodicLinks, numAffine, numCorrespondingNodes, *nodeTags = NULL, link, node;
848:   double         dbuf[16];

852:   GmshReadSize(gmsh, &numPeriodicLinks, 1);
853:   for (link = 0; link < numPeriodicLinks; ++link) {
854:     GmshReadInt(gmsh, info, 3);
855:     GmshReadSize(gmsh, &numAffine, 1);
856:     GmshReadDouble(gmsh, dbuf, numAffine);
857:     GmshReadSize(gmsh, &numCorrespondingNodes, 1);
858:     GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags);
859:     GmshReadSize(gmsh, nodeTags, numCorrespondingNodes*2);
860:     for (node = 0; node < numCorrespondingNodes; ++node) {
861:       PetscInt slaveNode  = nodeTags[node*2+0] - shift;
862:       PetscInt masterNode = nodeTags[node*2+1] - shift;
863:       slaveMap[slaveNode] = masterNode;
864:       PetscBTSet(bt, slaveNode);
865:       PetscBTSet(bt, masterNode);
866:     }
867:   }
868:   return(0);
869: }

871: /*
872: $MeshFormat // same as MSH version 2
873:   version(ASCII double; currently 4.1)
874:   file-type(ASCII int; 0 for ASCII mode, 1 for binary mode)
875:   data-size(ASCII int; sizeof(size_t))
876:   < int with value one; only in binary mode, to detect endianness >
877: $EndMeshFormat
878: */
879: static PetscErrorCode DMPlexCreateGmsh_ReadMeshFormat(GmshFile *gmsh)
880: {
881:   char           line[PETSC_MAX_PATH_LEN];
882:   int            snum, fileType, fileFormat, dataSize, checkEndian;
883:   float          version;

887:   GmshReadString(gmsh, line, 3);
888:   snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
889:   if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
890:   if (version < 2.2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at least 2.2", (double)version);
891:   if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f not supported", (double)version);
892:   if (version > 4.1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at most 4.1", (double)version);
893:   if (gmsh->binary && !fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Viewer is binary but Gmsh file is ASCII");
894:   if (!gmsh->binary && fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Viewer is ASCII but Gmsh file is binary");
895:   fileFormat = (int)(version*10); /* XXX Should use (int)roundf(version*10) ? */
896:   if (fileFormat <= 40 && dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
897:   if (fileFormat >= 41 && dataSize != sizeof(int) && dataSize != sizeof(PetscInt64)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
898:   gmsh->fileFormat = fileFormat;
899:   gmsh->dataSize = dataSize;
900:   gmsh->byteSwap = PETSC_FALSE;
901:   if (gmsh->binary) {
902:     GmshReadInt(gmsh, &checkEndian, 1);
903:     if (checkEndian != 1) {
904:       PetscByteSwap(&checkEndian, PETSC_ENUM, 1);
905:       if (checkEndian != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to detect endianness in Gmsh file header: %s", line);
906:       gmsh->byteSwap = PETSC_TRUE;
907:     }
908:   }
909:   return(0);
910: }

912: /*
913: PhysicalNames
914:   numPhysicalNames(ASCII int)
915:   dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max)
916:   ...
917: $EndPhysicalNames
918: */
919: static PetscErrorCode DMPlexCreateGmsh_ReadPhysicalNames(GmshFile *gmsh)
920: {
921:   char           line[PETSC_MAX_PATH_LEN], name[128+2], *p, *q;
922:   int            snum, numRegions, region, dim, tag;

926:   GmshReadString(gmsh, line, 1);
927:   snum = sscanf(line, "%d", &numRegions);
928:   if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
929:   for (region = 0; region < numRegions; ++region) {
930:     GmshReadString(gmsh, line, 2);
931:     snum = sscanf(line, "%d %d", &dim, &tag);
932:     if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
933:     GmshReadString(gmsh, line, -(PetscInt)sizeof(line));
934:     PetscStrchr(line, '"', &p);
935:     if (!p) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
936:     PetscStrrchr(line, '"', &q);
937:     if (q == p) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
938:     PetscStrncpy(name, p+1, (size_t)(q-p-1));
939:   }
940:   return(0);
941: }


944: PETSC_STATIC_INLINE DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt ctGmsh)
945: {
946:   switch (ctGmsh) {
947:     case 15:
948:       return DM_POLYTOPE_POINT;
949:     case 1:
950:     case 8:
951:       return DM_POLYTOPE_SEGMENT;
952:     case 2:
953:     case 9:
954:       return DM_POLYTOPE_TRIANGLE;
955:     case 3:
956:     case 10:
957:       return DM_POLYTOPE_QUADRILATERAL;
958:     case 4:
959:     case 11:
960:       return DM_POLYTOPE_TETRAHEDRON;
961:     case 5:
962:     case 12:
963:       return DM_POLYTOPE_HEXAHEDRON;
964:     case 6:
965:     case 13:
966:       return DM_POLYTOPE_TRI_PRISM;
967:     case 7:
968:     case 14:
969:       return DM_POLYTOPE_UNKNOWN; /* Pyramid */
970:     default:
971:       return DM_POLYTOPE_UNKNOWN;
972:   }
973: }

975: /*@C
976:   DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file

978: + comm        - The MPI communicator
979: . filename    - Name of the Gmsh file
980: - interpolate - Create faces and edges in the mesh

982:   Output Parameter:
983: . dm  - The DM object representing the mesh

985:   Level: beginner

987: .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
988: @*/
989: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
990: {
991:   PetscViewer     viewer;
992:   PetscMPIInt     rank;
993:   int             fileType;
994:   PetscViewerType vtype;
995:   PetscErrorCode  ierr;

998:   MPI_Comm_rank(comm, &rank);

1000:   /* Determine Gmsh file type (ASCII or binary) from file header */
1001:   if (!rank) {
1002:     GmshFile    gmsh_, *gmsh = &gmsh_;
1003:     char        line[PETSC_MAX_PATH_LEN];
1004:     int         snum;
1005:     float       version;

1007:     PetscArrayzero(gmsh,1);
1008:     PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer);
1009:     PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII);
1010:     PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ);
1011:     PetscViewerFileSetName(gmsh->viewer, filename);
1012:     /* Read only the first two lines of the Gmsh file */
1013:     GmshReadSection(gmsh, line);
1014:     GmshExpect(gmsh, "$MeshFormat", line);
1015:     GmshReadString(gmsh, line, 2);
1016:     snum = sscanf(line, "%f %d", &version, &fileType);
1017:     if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
1018:     if (version < 2.2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at least 2.2", (double)version);
1019:     if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f not supported", (double)version);
1020:     if (version > 4.1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at most 4.1", (double)version);
1021:     PetscViewerDestroy(&gmsh->viewer);
1022:   }
1023:   MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);
1024:   vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;

1026:   /* Create appropriate viewer and build plex */
1027:   PetscViewerCreate(comm, &viewer);
1028:   PetscViewerSetType(viewer, vtype);
1029:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);
1030:   PetscViewerFileSetName(viewer, filename);
1031:   DMPlexCreateGmsh(comm, viewer, interpolate, dm);
1032:   PetscViewerDestroy(&viewer);
1033:   return(0);
1034: }

1036: /*@
1037:   DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer

1039:   Collective

1041:   Input Parameters:
1042: + comm  - The MPI communicator
1043: . viewer - The Viewer associated with a Gmsh file
1044: - interpolate - Create faces and edges in the mesh

1046:   Output Parameter:
1047: . dm  - The DM object representing the mesh

1049:   Note: http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format

1051:   Level: beginner

1053: .seealso: DMPLEX, DMCreate()
1054: @*/
1055: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1056: {
1057:   PetscViewer    parentviewer = NULL;
1058:   double        *coordsIn = NULL;
1059:   GmshEntities  *entities = NULL;
1060:   GmshElement   *gmsh_elem = NULL;
1061:   PetscSection   coordSection;
1062:   Vec            coordinates;
1063:   PetscBT        periodicV = NULL, periodicC = NULL;
1064:   PetscScalar   *coords;
1065:   PetscInt       dim = 0, embedDim = -1, coordSize, c, v, d, cell, *periodicMap = NULL, *periodicMapI = NULL, *hybridMap = NULL;
1066:   PetscInt       numVertices = 0, numCells = 0, trueNumCells = 0;
1067:   int            i, shift = 1;
1068:   PetscMPIInt    rank;
1069:   PetscBool      binary, zerobase = PETSC_FALSE, usemarker = PETSC_FALSE;
1070:   PetscBool      enable_hybrid = interpolate, periodic = PETSC_TRUE;
1071:   PetscBool      hasTetra = PETSC_FALSE;

1075:   MPI_Comm_rank(comm, &rank);
1076:   PetscObjectOptionsBegin((PetscObject)viewer);
1077:   PetscOptionsHead(PetscOptionsObject,"DMPlex Gmsh options");
1078:   PetscOptionsBool("-dm_plex_gmsh_hybrid", "Generate hybrid cell bounds", "DMPlexCreateGmsh", enable_hybrid, &enable_hybrid, NULL);
1079:   PetscOptionsBool("-dm_plex_gmsh_periodic","Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL);
1080:   PetscOptionsBool("-dm_plex_gmsh_use_marker", "Generate marker label", "DMPlexCreateGmsh", usemarker, &usemarker, NULL);
1081:   PetscOptionsBool("-dm_plex_gmsh_zero_base", "Read Gmsh file with zero base indices", "DMPlexCreateGmsh", zerobase, &zerobase, NULL);
1082:   PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", embedDim, &embedDim, NULL,PETSC_DECIDE);
1083:   PetscOptionsTail();
1084:   PetscOptionsEnd();
1085:   if (zerobase) shift = 0;

1087:   DMCreate(comm, dm);
1088:   DMSetType(*dm, DMPLEX);
1089:   PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);

1091:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);

1093:   /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1094:   if (binary) {
1095:     parentviewer = viewer;
1096:     PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1097:   }

1099:   if (!rank) {
1100:     GmshFile  gmsh_, *gmsh = &gmsh_;
1101:     char      line[PETSC_MAX_PATH_LEN];
1102:     PetscBool match;

1104:     PetscArrayzero(gmsh,1);
1105:     gmsh->viewer = viewer;
1106:     gmsh->binary = binary;

1108:     /* Read mesh format */
1109:     GmshReadSection(gmsh, line);
1110:     GmshExpect(gmsh, "$MeshFormat", line);
1111:     DMPlexCreateGmsh_ReadMeshFormat(gmsh);
1112:     GmshReadEndSection(gmsh, "$EndMeshFormat", line);

1114:     /* OPTIONAL Read physical names */
1115:     GmshReadSection(gmsh, line);
1116:     GmshMatch(gmsh,"$PhysicalNames", line, &match);
1117:     if (match) {
1118:       DMPlexCreateGmsh_ReadPhysicalNames(gmsh);
1119:       GmshReadEndSection(gmsh, "$EndPhysicalNames", line);
1120:       /* Initial read for entity section */
1121:       GmshReadSection(gmsh, line);
1122:     }

1124:     /* Read entities */
1125:     if (gmsh->fileFormat >= 40) {
1126:       GmshExpect(gmsh, "$Entities", line);
1127:       switch (gmsh->fileFormat) {
1128:       case 41: DMPlexCreateGmsh_ReadEntities_v41(gmsh, &entities); break;
1129:       default: DMPlexCreateGmsh_ReadEntities_v40(gmsh, &entities); break;
1130:       }
1131:       GmshReadEndSection(gmsh, "$EndEntities", line);
1132:       /* Initial read for nodes section */
1133:       GmshReadSection(gmsh, line);
1134:     }

1136:     /* Read nodes */
1137:     GmshExpect(gmsh, "$Nodes", line);
1138:     switch (gmsh->fileFormat) {
1139:     case 41: DMPlexCreateGmsh_ReadNodes_v41(gmsh, shift, &numVertices, &coordsIn); break;
1140:     case 40: DMPlexCreateGmsh_ReadNodes_v40(gmsh, shift, &numVertices, &coordsIn); break;
1141:     default: DMPlexCreateGmsh_ReadNodes_v22(gmsh, shift, &numVertices, &coordsIn); break;
1142:     }
1143:     GmshReadEndSection(gmsh, "$EndNodes", line);

1145:     /* Read elements */
1146:     GmshReadSection(gmsh, line);
1147:     GmshExpect(gmsh, "$Elements", line);
1148:     switch (gmsh->fileFormat) {
1149:     case 41: DMPlexCreateGmsh_ReadElements_v41(gmsh, shift, entities, &numCells, &gmsh_elem); break;
1150:     case 40: DMPlexCreateGmsh_ReadElements_v40(gmsh, shift, entities, &numCells, &gmsh_elem); break;
1151:     default: DMPlexCreateGmsh_ReadElements_v22(gmsh, shift, &numCells, &gmsh_elem); break;
1152:     }
1153:     GmshReadEndSection(gmsh, "$EndElements", line);

1155:     /* OPTIONAL Read periodic section */
1156:     if (periodic) {
1157:       GmshReadSection(gmsh, line);
1158:       GmshMatch(gmsh, "$Periodic", line, &periodic);
1159:     }
1160:     if (periodic) {
1161:       PetscInt pVert, *periodicMapT, *aux;

1163:       PetscMalloc1(numVertices, &periodicMapT);
1164:       PetscBTCreate(numVertices, &periodicV);
1165:       for (i = 0; i < numVertices; i++) periodicMapT[i] = i;

1167:       GmshExpect(gmsh, "$Periodic", line);
1168:       switch (gmsh->fileFormat) {
1169:       case 41: DMPlexCreateGmsh_ReadPeriodic_v41(gmsh, shift, periodicMapT, periodicV); break;
1170:       default: DMPlexCreateGmsh_ReadPeriodic_v40(gmsh, shift, periodicMapT, periodicV); break;
1171:       }
1172:       GmshReadEndSection(gmsh, "$EndPeriodic", line);

1174:       /* we may have slaves of slaves */
1175:       for (i = 0; i < numVertices; i++) {
1176:         while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) {
1177:           periodicMapT[i] = periodicMapT[periodicMapT[i]];
1178:         }
1179:       }
1180:       /* periodicMap : from old to new numbering (periodic vertices excluded)
1181:          periodicMapI: from new to old numbering */
1182:       PetscMalloc1(numVertices, &periodicMap);
1183:       PetscMalloc1(numVertices, &periodicMapI);
1184:       PetscMalloc1(numVertices, &aux);
1185:       for (i = 0, pVert = 0; i < numVertices; i++) {
1186:         if (periodicMapT[i] != i) {
1187:           pVert++;
1188:         } else {
1189:           aux[i] = i - pVert;
1190:           periodicMapI[i - pVert] = i;
1191:         }
1192:       }
1193:       for (i = 0 ; i < numVertices; i++) {
1194:         periodicMap[i] = aux[periodicMapT[i]];
1195:       }
1196:       PetscFree(periodicMapT);
1197:       PetscFree(aux);
1198:       /* remove periodic vertices */
1199:       numVertices = numVertices - pVert;
1200:     }

1202:     GmshEntitiesDestroy(&entities);
1203:     PetscFree(gmsh->wbuf);
1204:     PetscFree(gmsh->sbuf);
1205:   }

1207:   if (parentviewer) {
1208:     PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1209:   }

1211:   if (!rank) {
1212:     PetscBool hybrid   = PETSC_FALSE;
1213:     PetscInt  cellType = -1;

1215:     for (trueNumCells = 0, c = 0; c < numCells; ++c) {
1216:       if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0; hybrid = PETSC_FALSE; cellType = -1;}
1217:       if (gmsh_elem[c].dim < dim) continue;
1218:       if (cellType == -1) cellType = gmsh_elem[c].cellType;
1219:       /* different cell type indicate an hybrid mesh in PLEX */
1220:       if (cellType != gmsh_elem[c].cellType) hybrid = PETSC_TRUE;
1221:       /* wedges always indicate an hybrid mesh in PLEX */
1222:       if (cellType == 6 || cellType == 13) hybrid = PETSC_TRUE;
1223:       if (cellType == 4 || cellType == 11) hasTetra = PETSC_TRUE;
1224:       trueNumCells++;
1225:     }
1226:     /* Renumber cells for hybrid grids */
1227:     if (hybrid && enable_hybrid) {
1228:       PetscInt hc1 = 0, hc2 = 0, *hybridCells1 = NULL, *hybridCells2 = NULL;
1229:       PetscInt cell, tn, *tp;
1230:       int      n1 = 0,n2 = 0;

1232:       PetscMalloc1(trueNumCells, &hybridCells1);
1233:       PetscMalloc1(trueNumCells, &hybridCells2);
1234:       for (cell = 0, c = 0; c < numCells; ++c) {
1235:         if (gmsh_elem[c].dim == dim) {
1236:           if (!n1) n1 = gmsh_elem[c].cellType;
1237:           else if (!n2 && n1 != gmsh_elem[c].cellType) n2 = gmsh_elem[c].cellType;

1239:           if      (gmsh_elem[c].cellType == n1) hybridCells1[hc1++] = cell;
1240:           else if (gmsh_elem[c].cellType == n2) hybridCells2[hc2++] = cell;
1241:           else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle more than 2 cell types");
1242:           cell++;
1243:         }
1244:       }

1246:       switch (n1) {
1247:       case 2: /* triangles */
1248:       case 9:
1249:         switch (n2) {
1250:         case 0: /* single type mesh */
1251:         case 3: /* quads */
1252:         case 10:
1253:           break;
1254:         default:
1255:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1256:         }
1257:         break;
1258:       case 3: /* quadrilateral */
1259:       case 10:
1260:         switch (n2) {
1261:         case 0: /* single type mesh */
1262:         case 2: /* swap since we list simplices first */
1263:         case 9:
1264:           tn  = hc1;
1265:           hc1 = hc2;
1266:           hc2 = tn;

1268:           tp           = hybridCells1;
1269:           hybridCells1 = hybridCells2;
1270:           hybridCells2 = tp;
1271:           break;
1272:         default:
1273:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1274:         }
1275:         break;
1276:       case 4: /* tetrahedra */
1277:       case 11:
1278:         switch (n2) {
1279:         case 0: /* single type mesh */
1280:         case 6: /* wedges */
1281:         case 13:
1282:           break;
1283:         default:
1284:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1285:         }
1286:         break;
1287:       case 5: /* hexahedra */
1288:       case 12:
1289:         switch (n2) {
1290:         case 0: /* single type mesh */
1291:         case 6: /* wedges */
1292:         case 13:
1293:           break;
1294:         default:
1295:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1296:         }
1297:         break;
1298:       case 6: /* wedge */
1299:       case 13:
1300:         switch (n2) {
1301:         case 0: /* single type mesh */
1302:         case 4: /* tetrahedra: swap since we list simplices first */
1303:         case 11:
1304:         case 5: /* hexahedra */
1305:         case 12:
1306:           tn  = hc1;
1307:           hc1 = hc2;
1308:           hc2 = tn;

1310:           tp           = hybridCells1;
1311:           hybridCells1 = hybridCells2;
1312:           hybridCells2 = tp;
1313:           break;
1314:         default:
1315:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1316:         }
1317:         break;
1318:       default:
1319:         SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1320:       }
1321:       PetscMalloc1(trueNumCells, &hybridMap);
1322:       for (cell = 0; cell < hc1; cell++) hybridMap[hybridCells1[cell]] = cell;
1323:       for (cell = 0; cell < hc2; cell++) hybridMap[hybridCells2[cell]] = cell + hc1;
1324:       PetscFree(hybridCells1);
1325:       PetscFree(hybridCells2);
1326:     }

1328:   }

1330:   /* Allocate the cell-vertex mesh */
1331:   /*   We do not want this label automatically computed, instead we compute it here */
1332:   DMCreateLabel(*dm, "celltype");
1333:   DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
1334:   for (cell = 0, c = 0; c < numCells; ++c) {
1335:     if (gmsh_elem[c].dim == dim) {
1336:       PetscInt       ucell = hybridMap ? hybridMap[cell] : cell;
1337:       DMPolytopeType ctype = DMPolytopeTypeFromGmsh(gmsh_elem[c].cellType);
1338:       if (hybridMap && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR;
1339:       DMPlexSetConeSize(*dm, ucell, gmsh_elem[c].numNodes);
1340:       DMPlexSetCellType(*dm, ucell, ctype);
1341:       cell++;
1342:     }
1343:   }
1344:   for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
1345:     DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
1346:   }
1347:   DMSetUp(*dm);

1349:   /* Add cell-vertex connections */
1350:   for (cell = 0, c = 0; c < numCells; ++c) {
1351:     if (gmsh_elem[c].dim == dim) {
1352:       PetscInt pcone[8], corner;
1353:       PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1354:       for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1355:         const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
1356:         pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells;
1357:       }
1358:       DMPlexReorderCell(*dm, ucell, pcone);
1359:       DMPlexSetCone(*dm, ucell, pcone);
1360:       cell++;
1361:     }
1362:   }

1364:   MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
1365:   DMSetDimension(*dm, dim);
1366:   DMPlexSymmetrize(*dm);
1367:   DMPlexStratify(*dm);
1368:   if (interpolate) {
1369:     DM idm;

1371:     DMPlexInterpolate(*dm, &idm);
1372:     DMDestroy(dm);
1373:     *dm  = idm;
1374:   }

1376:   if (usemarker && !interpolate && dim > 1) SETERRQ(comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation");
1377:   if (!rank && usemarker) {
1378:     PetscInt f, fStart, fEnd;

1380:     DMCreateLabel(*dm, "marker");
1381:     DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
1382:     for (f = fStart; f < fEnd; ++f) {
1383:       PetscInt suppSize;

1385:       DMPlexGetSupportSize(*dm, f, &suppSize);
1386:       if (suppSize == 1) {
1387:         PetscInt *cone = NULL, coneSize, p;

1389:         DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1390:         for (p = 0; p < coneSize; p += 2) {
1391:           DMSetLabelValue(*dm, "marker", cone[p], 1);
1392:         }
1393:         DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1394:       }
1395:     }
1396:   }

1398:   if (!rank) {
1399:     PetscInt vStart, vEnd;

1401:     DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
1402:     for (cell = 0, c = 0; c < numCells; ++c) {

1404:       /* Create face sets */
1405:       if (interpolate && gmsh_elem[c].dim == dim-1) {
1406:         const PetscInt *join;
1407:         PetscInt        joinSize, pcone[8], corner;
1408:         /* Find the relevant facet with vertex joins */
1409:         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1410:           const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
1411:           pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart;
1412:         }
1413:         DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, pcone, &joinSize, &join);
1414:         if (joinSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for GMsh element %d (Plex cell %D)", gmsh_elem[c].id, c);
1415:         DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
1416:         DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
1417:       }

1419:       /* Create cell sets */
1420:       if (gmsh_elem[c].dim == dim) {
1421:         if (gmsh_elem[c].numTags > 0) {
1422:           DMSetLabelValue(*dm, "Cell Sets", hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].tags[0]);
1423:         }
1424:         cell++;
1425:       }

1427:       /* Create vertex sets */
1428:       if (gmsh_elem[c].dim == 0) {
1429:         if (gmsh_elem[c].numTags > 0) {
1430:           const PetscInt cc = gmsh_elem[c].nodes[0] - shift;
1431:           const PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart;
1432:           DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);
1433:         }
1434:       }
1435:     }
1436:   }

1438:   /* Create coordinates */
1439:   if (embedDim < 0) embedDim = dim;
1440:   DMSetCoordinateDim(*dm, embedDim);
1441:   DMGetCoordinateSection(*dm, &coordSection);
1442:   PetscSectionSetNumFields(coordSection, 1);
1443:   PetscSectionSetFieldComponents(coordSection, 0, embedDim);
1444:   if (periodicMap) { /* we need to localize coordinates on cells */
1445:     PetscSectionSetChart(coordSection, 0, trueNumCells + numVertices);
1446:   } else {
1447:     PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);
1448:   }
1449:   for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
1450:     PetscSectionSetDof(coordSection, v, embedDim);
1451:     PetscSectionSetFieldDof(coordSection, v, 0, embedDim);
1452:   }
1453:   if (periodicMap) {
1454:     PetscBTCreate(trueNumCells, &periodicC);
1455:     for (cell = 0, c = 0; c < numCells; ++c) {
1456:       if (gmsh_elem[c].dim == dim) {
1457:         PetscInt  corner, pc = PETSC_FALSE;
1458:         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1459:           pc = (PetscBool)(pc || PetscBTLookup(periodicV, gmsh_elem[c].nodes[corner] - shift));
1460:         }
1461:         if (pc) {
1462:           PetscInt dof = gmsh_elem[c].numNodes*embedDim;
1463:           PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1464:           PetscBTSet(periodicC, ucell);
1465:           PetscSectionSetDof(coordSection, ucell, dof);
1466:           PetscSectionSetFieldDof(coordSection, ucell, 0, dof);
1467:         }
1468:         cell++;
1469:       }
1470:     }
1471:   }
1472:   PetscSectionSetUp(coordSection);
1473:   PetscSectionGetStorageSize(coordSection, &coordSize);
1474:   VecCreate(PETSC_COMM_SELF, &coordinates);
1475:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
1476:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
1477:   VecSetBlockSize(coordinates, embedDim);
1478:   VecSetType(coordinates, VECSTANDARD);
1479:   VecGetArray(coordinates, &coords);
1480:   if (periodicMap) {
1481:     PetscInt off;

1483:     for (cell = 0, c = 0; c < numCells; ++c) {
1484:       if (gmsh_elem[c].dim == dim) {
1485:         PetscInt pcone[8], corner;
1486:         PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1487:         if (PetscUnlikely(PetscBTLookup(periodicC, ucell))) {
1488:           for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1489:             pcone[corner] = gmsh_elem[c].nodes[corner] - shift;
1490:           }
1491:           DMPlexReorderCell(*dm, ucell, pcone);
1492:           PetscSectionGetOffset(coordSection, ucell, &off);
1493:           for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner)
1494:             for (v = pcone[corner], d = 0; d < embedDim; ++d)
1495:               coords[off++] = (PetscReal) coordsIn[v*3+d];
1496:         }
1497:         cell++;
1498:       }
1499:     }
1500:     for (v = 0; v < numVertices; ++v) {
1501:       PetscSectionGetOffset(coordSection, v + trueNumCells, &off);
1502:       for (d = 0; d < embedDim; ++d) {
1503:         coords[off+d] = (PetscReal) coordsIn[periodicMapI[v]*3+d];
1504:       }
1505:     }
1506:   } else {
1507:     for (v = 0; v < numVertices; ++v) {
1508:       for (d = 0; d < embedDim; ++d) {
1509:         coords[v*embedDim+d] = (PetscReal) coordsIn[v*3+d];
1510:       }
1511:     }
1512:   }
1513:   VecRestoreArray(coordinates, &coords);
1514:   DMSetCoordinatesLocal(*dm, coordinates);

1516:   periodic = periodicMap ? PETSC_TRUE : PETSC_FALSE;
1517:   MPI_Bcast(&periodic, 1, MPIU_BOOL, 0, comm);
1518:   DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);

1520:   PetscFree(coordsIn);
1521:   PetscFree(gmsh_elem);
1522:   PetscFree(hybridMap);
1523:   PetscFree(periodicMap);
1524:   PetscFree(periodicMapI);
1525:   PetscBTDestroy(&periodicV);
1526:   PetscBTDestroy(&periodicC);
1527:   VecDestroy(&coordinates);

1529:   PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
1530:   return(0);
1531: }