Actual source code: plexgmsh.c

petsc-3.11.4 2019-09-28
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 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: }