Actual source code: plexgmsh.c

petsc-master 2020-08-05
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:   PetscInt     nodeTagMin, nodeTagMax;
 16:   PetscInt    *nodeMap;
 17: } GmshFile;

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

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

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

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

 50: static PetscErrorCode GmshRead(GmshFile *gmsh, void *buf, PetscInt count, PetscDataType dtype)
 51: {
 54:   PetscViewerRead(gmsh->viewer, buf, count, NULL, dtype);
 55:   if (gmsh->byteSwap) {PetscByteSwap(buf, dtype, count);}
 56:   return(0);
 57: }

 59: static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscInt count)
 60: {
 63:   PetscViewerRead(gmsh->viewer, buf, count, NULL, PETSC_STRING);
 64:   return(0);
 65: }

 67: static PetscErrorCode GmshMatch(PETSC_UNUSED GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN], PetscBool *match)
 68: {
 71:   PetscStrcmp(line, Section, match);
 72:   return(0);
 73: }

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

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

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

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

105: static PetscErrorCode GmshReadEndSection(GmshFile *gmsh, const char EndSection[], char line[PETSC_MAX_PATH_LEN])
106: {
109:   GmshReadString(gmsh, line, 1);
110:   GmshExpect(gmsh, EndSection, line);
111:   return(0);
112: }

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

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

142: static PetscErrorCode GmshReadInt(GmshFile *gmsh, int *buf, PetscInt count)
143: {
146:   GmshRead(gmsh, buf, count, PETSC_ENUM);
147:   return(0);
148: }

150: static PetscErrorCode GmshReadDouble(GmshFile *gmsh, double *buf, PetscInt count)
151: {
154:   GmshRead(gmsh, buf, count, PETSC_DOUBLE);
155:   return(0);
156: }

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

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

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

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

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

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

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

207: static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities)
208: {
209:   PetscInt       dim;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

723: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
724: $Nodes
725:   numEntityBlocks(size_t) numNodes(size_t)
726:     minNodeTag(size_t) maxNodeTag(size_t)
727:   entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesBlock(size_t)
728:     nodeTag(size_t)
729:     ...
730:     x(double) y(double) z(double)
731:        < u(double; if parametric and entityDim = 1 or entityDim = 2) >
732:        < v(double; if parametric and entityDim = 2) >
733:     ...
734:   ...
735: $EndNodes

737: We need to allow holes in the node numbering, so we must have a map between the file numbering and a contiguous numbering.
738: */
739: static PetscErrorCode DMPlexCreateGmsh_ReadNodes_v41(GmshFile *gmsh, int shift, PetscInt *numVertices, double **gmsh_nodes)
740: {
741:   int            info[3];
742:   PetscInt       sizes[4], numEntityBlocks, numNodes, numNodesBlock = 0, block, *nodeTag = NULL, node, i;
743:   double         *coordinates;

747:   GmshReadSize(gmsh, sizes, 4);
748:   numEntityBlocks = sizes[0]; numNodes = sizes[1]; gmsh->nodeTagMin = sizes[2]; gmsh->nodeTagMax = sizes[3]+1;
749:   PetscMalloc1(numNodes*3, &coordinates);
750:   PetscMalloc1(gmsh->nodeTagMax-gmsh->nodeTagMin, &gmsh->nodeMap);
751:   for (i = 0; i < gmsh->nodeTagMax-gmsh->nodeTagMin; ++i) gmsh->nodeMap[i] = -1;
752:   *numVertices = numNodes;
753:   *gmsh_nodes = coordinates;
754:   for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) {
755:     GmshReadInt(gmsh, info, 3);
756:     GmshReadSize(gmsh, &numNodesBlock, 1);
757:     GmshBufferGet(gmsh, numNodesBlock, sizeof(PetscInt), &nodeTag);
758:     GmshReadSize(gmsh, nodeTag, numNodesBlock);
759:     for (i = 0; i < numNodesBlock; ++i) {
760:       const PetscInt off = nodeTag[i] - gmsh->nodeTagMin;

762:       if (gmsh->nodeMap[off] != -1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Repeated node tag %D", nodeTag[i]);
763:       gmsh->nodeMap[off] = node+i+shift;
764:     }
765:     if (info[2] != 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parametric coordinates not supported");
766:     GmshReadDouble(gmsh, coordinates+node*3, numNodesBlock*3);
767:   }
768:   return(0);
769: }

771: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
772: $Elements
773:   numEntityBlocks(size_t) numElements(size_t)
774:     minElementTag(size_t) maxElementTag(size_t)
775:   entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t)
776:     elementTag(size_t) nodeTag(size_t) ...
777:     ...
778:   ...
779: $EndElements
780: */
781: static PetscErrorCode DMPlexCreateGmsh_ReadElements_v41(GmshFile *gmsh, PETSC_UNUSED int shift, GmshEntities *entities, PetscInt *numCells, GmshElement **gmsh_elems)
782: {
783:   int            info[3], eid, dim, cellType, *tags;
784:   PetscInt       nodeTagMin = gmsh->nodeTagMin, *nodeMap = gmsh->nodeMap, numNodes;
785:   PetscInt       sizes[4], *ibuf = NULL, numEntityBlocks, numElements, numBlockElements, numTags, block, elem, c, p;
786:   GmshEntity     *entity = NULL;
787:   GmshElement    *elements;

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

845: /*
846: $Periodic
847:   numPeriodicLinks(size_t)
848:  entityDim(int) entityTag(int) entityTagMaster(int)
849:   numAffine(size_t) value(double) ...
850:   numCorrespondingNodes(size_t)
851:     nodeTag(size_t) nodeTagMaster(size_t)
852:     ...
853:   ...
854: $EndPeriodic
855: */
856: static PetscErrorCode DMPlexCreateGmsh_ReadPeriodic_v41(GmshFile *gmsh, int shift, PetscInt slaveMap[], PetscBT bt)
857: {
858:   int            info[3];
859:   PetscInt       numPeriodicLinks, numAffine, numCorrespondingNodes, *nodeTags = NULL, link, node;
860:   double         dbuf[16];

864:   GmshReadSize(gmsh, &numPeriodicLinks, 1);
865:   for (link = 0; link < numPeriodicLinks; ++link) {
866:     GmshReadInt(gmsh, info, 3);
867:     GmshReadSize(gmsh, &numAffine, 1);
868:     GmshReadDouble(gmsh, dbuf, numAffine);
869:     GmshReadSize(gmsh, &numCorrespondingNodes, 1);
870:     GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags);
871:     GmshReadSize(gmsh, nodeTags, numCorrespondingNodes*2);
872:     for (node = 0; node < numCorrespondingNodes; ++node) {
873:       PetscInt slaveNode  = nodeTags[node*2+0] - shift;
874:       PetscInt masterNode = nodeTags[node*2+1] - shift;
875:       slaveMap[slaveNode] = masterNode;
876:       PetscBTSet(bt, slaveNode);
877:       PetscBTSet(bt, masterNode);
878:     }
879:   }
880:   return(0);
881: }

883: /*
884: $MeshFormat // same as MSH version 2
885:   version(ASCII double; currently 4.1)
886:   file-type(ASCII int; 0 for ASCII mode, 1 for binary mode)
887:   data-size(ASCII int; sizeof(size_t))
888:   < int with value one; only in binary mode, to detect endianness >
889: $EndMeshFormat
890: */
891: static PetscErrorCode DMPlexCreateGmsh_ReadMeshFormat(GmshFile *gmsh)
892: {
893:   char           line[PETSC_MAX_PATH_LEN];
894:   int            snum, fileType, fileFormat, dataSize, checkEndian;
895:   float          version;

899:   GmshReadString(gmsh, line, 3);
900:   snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
901:   if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
902:   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);
903:   if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f not supported", (double)version);
904:   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);
905:   if (gmsh->binary && !fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Viewer is binary but Gmsh file is ASCII");
906:   if (!gmsh->binary && fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Viewer is ASCII but Gmsh file is binary");
907:   fileFormat = (int)(version*10); /* XXX Should use (int)roundf(version*10) ? */
908:   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);
909:   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);
910:   gmsh->fileFormat = fileFormat;
911:   gmsh->dataSize = dataSize;
912:   gmsh->byteSwap = PETSC_FALSE;
913:   if (gmsh->binary) {
914:     GmshReadInt(gmsh, &checkEndian, 1);
915:     if (checkEndian != 1) {
916:       PetscByteSwap(&checkEndian, PETSC_ENUM, 1);
917:       if (checkEndian != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to detect endianness in Gmsh file header: %s", line);
918:       gmsh->byteSwap = PETSC_TRUE;
919:     }
920:   }
921:   return(0);
922: }

924: /*
925: PhysicalNames
926:   numPhysicalNames(ASCII int)
927:   dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max)
928:   ...
929: $EndPhysicalNames
930: */
931: static PetscErrorCode DMPlexCreateGmsh_ReadPhysicalNames(GmshFile *gmsh)
932: {
933:   char           line[PETSC_MAX_PATH_LEN], name[128+2], *p, *q;
934:   int            snum, numRegions, region, dim, tag;

938:   GmshReadString(gmsh, line, 1);
939:   snum = sscanf(line, "%d", &numRegions);
940:   if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
941:   for (region = 0; region < numRegions; ++region) {
942:     GmshReadString(gmsh, line, 2);
943:     snum = sscanf(line, "%d %d", &dim, &tag);
944:     if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
945:     GmshReadString(gmsh, line, -(PetscInt)sizeof(line));
946:     PetscStrchr(line, '"', &p);
947:     if (!p) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
948:     PetscStrrchr(line, '"', &q);
949:     if (q == p) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
950:     PetscStrncpy(name, p+1, (size_t)(q-p-1));
951:   }
952:   return(0);
953: }


956: PETSC_STATIC_INLINE DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt ctGmsh)
957: {
958:   switch (ctGmsh) {
959:     case 15:
960:       return DM_POLYTOPE_POINT;
961:     case 1:
962:     case 8:
963:       return DM_POLYTOPE_SEGMENT;
964:     case 2:
965:     case 9:
966:       return DM_POLYTOPE_TRIANGLE;
967:     case 3:
968:     case 10:
969:       return DM_POLYTOPE_QUADRILATERAL;
970:     case 4:
971:     case 11:
972:       return DM_POLYTOPE_TETRAHEDRON;
973:     case 5:
974:     case 12:
975:       return DM_POLYTOPE_HEXAHEDRON;
976:     case 6:
977:     case 13:
978:       return DM_POLYTOPE_TRI_PRISM;
979:     case 7:
980:     case 14:
981:       return DM_POLYTOPE_UNKNOWN; /* Pyramid */
982:     default:
983:       return DM_POLYTOPE_UNKNOWN;
984:   }
985: }

987: /*@C
988:   DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file

990: + comm        - The MPI communicator
991: . filename    - Name of the Gmsh file
992: - interpolate - Create faces and edges in the mesh

994:   Output Parameter:
995: . dm  - The DM object representing the mesh

997:   Level: beginner

999: .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
1000: @*/
1001: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1002: {
1003:   PetscViewer     viewer;
1004:   PetscMPIInt     rank;
1005:   int             fileType;
1006:   PetscViewerType vtype;
1007:   PetscErrorCode  ierr;

1010:   MPI_Comm_rank(comm, &rank);

1012:   /* Determine Gmsh file type (ASCII or binary) from file header */
1013:   if (!rank) {
1014:     GmshFile    gmsh_, *gmsh = &gmsh_;
1015:     char        line[PETSC_MAX_PATH_LEN];
1016:     int         snum;
1017:     float       version;

1019:     PetscArrayzero(gmsh,1);
1020:     PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer);
1021:     PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII);
1022:     PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ);
1023:     PetscViewerFileSetName(gmsh->viewer, filename);
1024:     /* Read only the first two lines of the Gmsh file */
1025:     GmshReadSection(gmsh, line);
1026:     GmshExpect(gmsh, "$MeshFormat", line);
1027:     GmshReadString(gmsh, line, 2);
1028:     snum = sscanf(line, "%f %d", &version, &fileType);
1029:     if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
1030:     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);
1031:     if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f not supported", (double)version);
1032:     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);
1033:     PetscViewerDestroy(&gmsh->viewer);
1034:   }
1035:   MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);
1036:   vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;

1038:   /* Create appropriate viewer and build plex */
1039:   PetscViewerCreate(comm, &viewer);
1040:   PetscViewerSetType(viewer, vtype);
1041:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);
1042:   PetscViewerFileSetName(viewer, filename);
1043:   DMPlexCreateGmsh(comm, viewer, interpolate, dm);
1044:   PetscViewerDestroy(&viewer);
1045:   return(0);
1046: }

1048: /*@
1049:   DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer

1051:   Collective

1053:   Input Parameters:
1054: + comm  - The MPI communicator
1055: . viewer - The Viewer associated with a Gmsh file
1056: - interpolate - Create faces and edges in the mesh

1058:   Output Parameter:
1059: . dm  - The DM object representing the mesh

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

1063:   Level: beginner

1065: .seealso: DMPLEX, DMCreate()
1066: @*/
1067: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1068: {
1069:   PetscViewer    parentviewer = NULL;
1070:   double        *coordsIn = NULL;
1071:   GmshEntities  *entities = NULL;
1072:   GmshElement   *gmsh_elem = NULL;
1073:   PetscSection   coordSection;
1074:   Vec            coordinates;
1075:   PetscBT        periodicV = NULL, periodicC = NULL;
1076:   PetscScalar   *coords;
1077:   PetscInt       dim = 0, embedDim = -1, coordSize, c, v, d, cell, *periodicMap = NULL, *periodicMapI = NULL, *hybridMap = NULL;
1078:   PetscInt       numVertices = 0, numCells = 0, trueNumCells = 0;
1079:   int            i, shift = 1;
1080:   PetscMPIInt    rank;
1081:   PetscBool      binary, zerobase = PETSC_FALSE, usemarker = PETSC_FALSE;
1082:   PetscBool      enable_hybrid = interpolate, periodic = PETSC_TRUE;
1083:   PetscBool      hasTetra = PETSC_FALSE;

1087:   MPI_Comm_rank(comm, &rank);
1088:   PetscObjectOptionsBegin((PetscObject)viewer);
1089:   PetscOptionsHead(PetscOptionsObject,"DMPlex Gmsh options");
1090:   PetscOptionsBool("-dm_plex_gmsh_hybrid", "Generate hybrid cell bounds", "DMPlexCreateGmsh", enable_hybrid, &enable_hybrid, NULL);
1091:   PetscOptionsBool("-dm_plex_gmsh_periodic","Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL);
1092:   PetscOptionsBool("-dm_plex_gmsh_use_marker", "Generate marker label", "DMPlexCreateGmsh", usemarker, &usemarker, NULL);
1093:   PetscOptionsBool("-dm_plex_gmsh_zero_base", "Read Gmsh file with zero base indices", "DMPlexCreateGmsh", zerobase, &zerobase, NULL);
1094:   PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", embedDim, &embedDim, NULL,PETSC_DECIDE);
1095:   PetscOptionsTail();
1096:   PetscOptionsEnd();
1097:   if (zerobase) shift = 0;

1099:   DMCreate(comm, dm);
1100:   DMSetType(*dm, DMPLEX);
1101:   PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);

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

1105:   /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1106:   if (binary) {
1107:     parentviewer = viewer;
1108:     PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1109:   }

1111:   if (!rank) {
1112:     GmshFile  gmsh_, *gmsh = &gmsh_;
1113:     char      line[PETSC_MAX_PATH_LEN];
1114:     PetscBool match;

1116:     PetscArrayzero(gmsh,1);
1117:     gmsh->viewer = viewer;
1118:     gmsh->binary = binary;

1120:     /* Read mesh format */
1121:     GmshReadSection(gmsh, line);
1122:     GmshExpect(gmsh, "$MeshFormat", line);
1123:     DMPlexCreateGmsh_ReadMeshFormat(gmsh);
1124:     GmshReadEndSection(gmsh, "$EndMeshFormat", line);

1126:     /* OPTIONAL Read physical names */
1127:     GmshReadSection(gmsh, line);
1128:     GmshMatch(gmsh,"$PhysicalNames", line, &match);
1129:     if (match) {
1130:       DMPlexCreateGmsh_ReadPhysicalNames(gmsh);
1131:       GmshReadEndSection(gmsh, "$EndPhysicalNames", line);
1132:       /* Initial read for entity section */
1133:       GmshReadSection(gmsh, line);
1134:     }

1136:     /* Read entities */
1137:     if (gmsh->fileFormat >= 40) {
1138:       GmshExpect(gmsh, "$Entities", line);
1139:       switch (gmsh->fileFormat) {
1140:       case 41: DMPlexCreateGmsh_ReadEntities_v41(gmsh, &entities); break;
1141:       default: DMPlexCreateGmsh_ReadEntities_v40(gmsh, &entities); break;
1142:       }
1143:       GmshReadEndSection(gmsh, "$EndEntities", line);
1144:       /* Initial read for nodes section */
1145:       GmshReadSection(gmsh, line);
1146:     }

1148:     /* Read nodes */
1149:     GmshExpect(gmsh, "$Nodes", line);
1150:     switch (gmsh->fileFormat) {
1151:     case 41: DMPlexCreateGmsh_ReadNodes_v41(gmsh, shift, &numVertices, &coordsIn); break;
1152:     case 40: DMPlexCreateGmsh_ReadNodes_v40(gmsh, shift, &numVertices, &coordsIn); break;
1153:     default: DMPlexCreateGmsh_ReadNodes_v22(gmsh, shift, &numVertices, &coordsIn); break;
1154:     }
1155:     GmshReadEndSection(gmsh, "$EndNodes", line);

1157:     /* Read elements */
1158:     GmshReadSection(gmsh, line);
1159:     GmshExpect(gmsh, "$Elements", line);
1160:     switch (gmsh->fileFormat) {
1161:     case 41: DMPlexCreateGmsh_ReadElements_v41(gmsh, shift, entities, &numCells, &gmsh_elem); break;
1162:     case 40: DMPlexCreateGmsh_ReadElements_v40(gmsh, shift, entities, &numCells, &gmsh_elem); break;
1163:     default: DMPlexCreateGmsh_ReadElements_v22(gmsh, shift, &numCells, &gmsh_elem); break;
1164:     }
1165:     GmshReadEndSection(gmsh, "$EndElements", line);

1167:     /* OPTIONAL Read periodic section */
1168:     if (periodic) {
1169:       GmshReadSection(gmsh, line);
1170:       GmshMatch(gmsh, "$Periodic", line, &periodic);
1171:     }
1172:     if (periodic) {
1173:       PetscInt pVert, *periodicMapT, *aux;

1175:       PetscMalloc1(numVertices, &periodicMapT);
1176:       PetscBTCreate(numVertices, &periodicV);
1177:       for (i = 0; i < numVertices; i++) periodicMapT[i] = i;

1179:       GmshExpect(gmsh, "$Periodic", line);
1180:       switch (gmsh->fileFormat) {
1181:       case 41: DMPlexCreateGmsh_ReadPeriodic_v41(gmsh, shift, periodicMapT, periodicV); break;
1182:       default: DMPlexCreateGmsh_ReadPeriodic_v40(gmsh, shift, periodicMapT, periodicV); break;
1183:       }
1184:       GmshReadEndSection(gmsh, "$EndPeriodic", line);

1186:       /* we may have slaves of slaves */
1187:       for (i = 0; i < numVertices; i++) {
1188:         while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) {
1189:           periodicMapT[i] = periodicMapT[periodicMapT[i]];
1190:         }
1191:       }
1192:       /* periodicMap : from old to new numbering (periodic vertices excluded)
1193:          periodicMapI: from new to old numbering */
1194:       PetscMalloc1(numVertices, &periodicMap);
1195:       PetscMalloc1(numVertices, &periodicMapI);
1196:       PetscMalloc1(numVertices, &aux);
1197:       for (i = 0, pVert = 0; i < numVertices; i++) {
1198:         if (periodicMapT[i] != i) {
1199:           pVert++;
1200:         } else {
1201:           aux[i] = i - pVert;
1202:           periodicMapI[i - pVert] = i;
1203:         }
1204:       }
1205:       for (i = 0 ; i < numVertices; i++) {
1206:         periodicMap[i] = aux[periodicMapT[i]];
1207:       }
1208:       PetscFree(periodicMapT);
1209:       PetscFree(aux);
1210:       /* remove periodic vertices */
1211:       numVertices = numVertices - pVert;
1212:     }

1214:     GmshEntitiesDestroy(&entities);
1215:     PetscFree(gmsh->wbuf);
1216:     PetscFree(gmsh->sbuf);
1217:     PetscFree(gmsh->nodeMap);
1218:   }

1220:   if (parentviewer) {
1221:     PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1222:   }

1224:   if (!rank) {
1225:     PetscBool hybrid   = PETSC_FALSE;
1226:     PetscInt  cellType = -1;

1228:     for (trueNumCells = 0, c = 0; c < numCells; ++c) {
1229:       if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0; hybrid = PETSC_FALSE; cellType = -1;}
1230:       if (gmsh_elem[c].dim < dim) continue;
1231:       if (cellType == -1) cellType = gmsh_elem[c].cellType;
1232:       /* different cell type indicate an hybrid mesh in PLEX */
1233:       if (cellType != gmsh_elem[c].cellType) hybrid = PETSC_TRUE;
1234:       /* wedges always indicate an hybrid mesh in PLEX */
1235:       if (cellType == 6 || cellType == 13) hybrid = PETSC_TRUE;
1236:       if (cellType == 4 || cellType == 11) hasTetra = PETSC_TRUE;
1237:       trueNumCells++;
1238:     }
1239:     /* Renumber cells for hybrid grids */
1240:     if (hybrid && enable_hybrid) {
1241:       PetscInt hc1 = 0, hc2 = 0, *hybridCells1 = NULL, *hybridCells2 = NULL;
1242:       PetscInt cell, tn, *tp;
1243:       int      n1 = 0,n2 = 0;

1245:       PetscMalloc1(trueNumCells, &hybridCells1);
1246:       PetscMalloc1(trueNumCells, &hybridCells2);
1247:       for (cell = 0, c = 0; c < numCells; ++c) {
1248:         if (gmsh_elem[c].dim == dim) {
1249:           if (!n1) n1 = gmsh_elem[c].cellType;
1250:           else if (!n2 && n1 != gmsh_elem[c].cellType) n2 = gmsh_elem[c].cellType;

1252:           if      (gmsh_elem[c].cellType == n1) hybridCells1[hc1++] = cell;
1253:           else if (gmsh_elem[c].cellType == n2) hybridCells2[hc2++] = cell;
1254:           else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle more than 2 cell types");
1255:           cell++;
1256:         }
1257:       }

1259:       switch (n1) {
1260:       case 2: /* triangles */
1261:       case 9:
1262:         switch (n2) {
1263:         case 0: /* single type mesh */
1264:         case 3: /* quads */
1265:         case 10:
1266:           break;
1267:         default:
1268:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1269:         }
1270:         break;
1271:       case 3: /* quadrilateral */
1272:       case 10:
1273:         switch (n2) {
1274:         case 0: /* single type mesh */
1275:         case 2: /* swap since we list simplices first */
1276:         case 9:
1277:           tn  = hc1;
1278:           hc1 = hc2;
1279:           hc2 = tn;

1281:           tp           = hybridCells1;
1282:           hybridCells1 = hybridCells2;
1283:           hybridCells2 = tp;
1284:           break;
1285:         default:
1286:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1287:         }
1288:         break;
1289:       case 4: /* tetrahedra */
1290:       case 11:
1291:         switch (n2) {
1292:         case 0: /* single type mesh */
1293:         case 6: /* wedges */
1294:         case 13:
1295:           break;
1296:         default:
1297:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1298:         }
1299:         break;
1300:       case 5: /* hexahedra */
1301:       case 12:
1302:         switch (n2) {
1303:         case 0: /* single type mesh */
1304:         case 6: /* wedges */
1305:         case 13:
1306:           break;
1307:         default:
1308:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1309:         }
1310:         break;
1311:       case 6: /* wedge */
1312:       case 13:
1313:         switch (n2) {
1314:         case 0: /* single type mesh */
1315:         case 4: /* tetrahedra: swap since we list simplices first */
1316:         case 11:
1317:         case 5: /* hexahedra */
1318:         case 12:
1319:           tn  = hc1;
1320:           hc1 = hc2;
1321:           hc2 = tn;

1323:           tp           = hybridCells1;
1324:           hybridCells1 = hybridCells2;
1325:           hybridCells2 = tp;
1326:           break;
1327:         default:
1328:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1329:         }
1330:         break;
1331:       default:
1332:         SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1333:       }
1334:       PetscMalloc1(trueNumCells, &hybridMap);
1335:       for (cell = 0; cell < hc1; cell++) hybridMap[hybridCells1[cell]] = cell;
1336:       for (cell = 0; cell < hc2; cell++) hybridMap[hybridCells2[cell]] = cell + hc1;
1337:       PetscFree(hybridCells1);
1338:       PetscFree(hybridCells2);
1339:     }

1341:   }

1343:   /* Allocate the cell-vertex mesh */
1344:   /*   We do not want this label automatically computed, instead we compute it here */
1345:   DMCreateLabel(*dm, "celltype");
1346:   DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
1347:   for (cell = 0, c = 0; c < numCells; ++c) {
1348:     if (gmsh_elem[c].dim == dim) {
1349:       PetscInt       ucell = hybridMap ? hybridMap[cell] : cell;
1350:       DMPolytopeType ctype = DMPolytopeTypeFromGmsh(gmsh_elem[c].cellType);
1351:       if (hybridMap && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR;
1352:       DMPlexSetConeSize(*dm, ucell, gmsh_elem[c].numNodes);
1353:       DMPlexSetCellType(*dm, ucell, ctype);
1354:       cell++;
1355:     }
1356:   }
1357:   for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
1358:     DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
1359:   }
1360:   DMSetUp(*dm);

1362:   /* Add cell-vertex connections */
1363:   for (cell = 0, c = 0; c < numCells; ++c) {
1364:     if (gmsh_elem[c].dim == dim) {
1365:       PetscInt pcone[8], corner;
1366:       PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1367:       for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1368:         const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
1369:         pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells;
1370:       }
1371:       DMPlexReorderCell(*dm, ucell, pcone);
1372:       DMPlexSetCone(*dm, ucell, pcone);
1373:       cell++;
1374:     }
1375:   }

1377:   MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
1378:   DMSetDimension(*dm, dim);
1379:   DMPlexSymmetrize(*dm);
1380:   DMPlexStratify(*dm);
1381:   if (interpolate) {
1382:     DM idm;

1384:     DMPlexInterpolate(*dm, &idm);
1385:     DMDestroy(dm);
1386:     *dm  = idm;
1387:   }

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

1393:     DMCreateLabel(*dm, "marker");
1394:     DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
1395:     for (f = fStart; f < fEnd; ++f) {
1396:       PetscInt suppSize;

1398:       DMPlexGetSupportSize(*dm, f, &suppSize);
1399:       if (suppSize == 1) {
1400:         PetscInt *cone = NULL, coneSize, p;

1402:         DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1403:         for (p = 0; p < coneSize; p += 2) {
1404:           DMSetLabelValue(*dm, "marker", cone[p], 1);
1405:         }
1406:         DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1407:       }
1408:     }
1409:   }

1411:   if (!rank) {
1412:     PetscInt vStart, vEnd;

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

1417:       /* Create face sets */
1418:       if (interpolate && gmsh_elem[c].dim == dim-1) {
1419:         const PetscInt *join;
1420:         PetscInt        joinSize, pcone[8], corner;
1421:         /* Find the relevant facet with vertex joins */
1422:         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1423:           const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
1424:           pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart;
1425:         }
1426:         DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, pcone, &joinSize, &join);
1427:         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);
1428:         DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
1429:         DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
1430:       }

1432:       /* Create cell sets */
1433:       if (gmsh_elem[c].dim == dim) {
1434:         if (gmsh_elem[c].numTags > 0) {
1435:           DMSetLabelValue(*dm, "Cell Sets", hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].tags[0]);
1436:         }
1437:         cell++;
1438:       }

1440:       /* Create vertex sets */
1441:       if (gmsh_elem[c].dim == 0) {
1442:         if (gmsh_elem[c].numTags > 0) {
1443:           const PetscInt cc = gmsh_elem[c].nodes[0] - shift;
1444:           const PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart;
1445:           DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);
1446:         }
1447:       }
1448:     }
1449:   }

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

1496:     for (cell = 0, c = 0; c < numCells; ++c) {
1497:       if (gmsh_elem[c].dim == dim) {
1498:         PetscInt pcone[8], corner;
1499:         PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1500:         if (PetscUnlikely(PetscBTLookup(periodicC, ucell))) {
1501:           for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1502:             pcone[corner] = gmsh_elem[c].nodes[corner] - shift;
1503:           }
1504:           DMPlexReorderCell(*dm, ucell, pcone);
1505:           PetscSectionGetOffset(coordSection, ucell, &off);
1506:           for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner)
1507:             for (v = pcone[corner], d = 0; d < embedDim; ++d)
1508:               coords[off++] = (PetscReal) coordsIn[v*3+d];
1509:         }
1510:         cell++;
1511:       }
1512:     }
1513:     for (v = 0; v < numVertices; ++v) {
1514:       PetscSectionGetOffset(coordSection, v + trueNumCells, &off);
1515:       for (d = 0; d < embedDim; ++d) {
1516:         coords[off+d] = (PetscReal) coordsIn[periodicMapI[v]*3+d];
1517:       }
1518:     }
1519:   } else {
1520:     for (v = 0; v < numVertices; ++v) {
1521:       for (d = 0; d < embedDim; ++d) {
1522:         coords[v*embedDim+d] = (PetscReal) coordsIn[v*3+d];
1523:       }
1524:     }
1525:   }
1526:   VecRestoreArray(coordinates, &coords);
1527:   DMSetCoordinatesLocal(*dm, coordinates);

1529:   periodic = periodicMap ? PETSC_TRUE : PETSC_FALSE;
1530:   MPI_Bcast(&periodic, 1, MPIU_BOOL, 0, comm);
1531:   DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);

1533:   PetscFree(coordsIn);
1534:   PetscFree(gmsh_elem);
1535:   PetscFree(hybridMap);
1536:   PetscFree(periodicMap);
1537:   PetscFree(periodicMapI);
1538:   PetscBTDestroy(&periodicV);
1539:   PetscBTDestroy(&periodicC);
1540:   VecDestroy(&coordinates);

1542:   PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
1543:   return(0);
1544: }