Actual source code: plexgmsh.c

petsc-3.12.0 2019-09-29
Report Typos and Errors
  1: #define PETSCDM_DLL
  2:  #include <petsc/private/dmpleximpl.h>
  3:  #include <petsc/private/hashmapi.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

943: /*@C
944:   DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file

946: + comm        - The MPI communicator
947: . filename    - Name of the Gmsh file
948: - interpolate - Create faces and edges in the mesh

950:   Output Parameter:
951: . dm  - The DM object representing the mesh

953:   Level: beginner

955: .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
956: @*/
957: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
958: {
959:   PetscViewer     viewer;
960:   PetscMPIInt     rank;
961:   int             fileType;
962:   PetscViewerType vtype;
963:   PetscErrorCode  ierr;

966:   MPI_Comm_rank(comm, &rank);

968:   /* Determine Gmsh file type (ASCII or binary) from file header */
969:   if (!rank) {
970:     GmshFile    gmsh_, *gmsh = &gmsh_;
971:     char        line[PETSC_MAX_PATH_LEN];
972:     int         snum;
973:     float       version;

975:     PetscArrayzero(gmsh,1);
976:     PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer);
977:     PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII);
978:     PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ);
979:     PetscViewerFileSetName(gmsh->viewer, filename);
980:     /* Read only the first two lines of the Gmsh file */
981:     GmshReadSection(gmsh, line);
982:     GmshExpect(gmsh, "$MeshFormat", line);
983:     GmshReadString(gmsh, line, 2);
984:     snum = sscanf(line, "%f %d", &version, &fileType);
985:     if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
986:     if (version < 2.2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at least 2.2", (double)version);
987:     if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f not supported", (double)version);
988:     if (version > 4.1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at most 4.1", (double)version);
989:     PetscViewerDestroy(&gmsh->viewer);
990:   }
991:   MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);
992:   vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;

994:   /* Create appropriate viewer and build plex */
995:   PetscViewerCreate(comm, &viewer);
996:   PetscViewerSetType(viewer, vtype);
997:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);
998:   PetscViewerFileSetName(viewer, filename);
999:   DMPlexCreateGmsh(comm, viewer, interpolate, dm);
1000:   PetscViewerDestroy(&viewer);
1001:   return(0);
1002: }

1004: /*@
1005:   DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer

1007:   Collective

1009:   Input Parameters:
1010: + comm  - The MPI communicator
1011: . viewer - The Viewer associated with a Gmsh file
1012: - interpolate - Create faces and edges in the mesh

1014:   Output Parameter:
1015: . dm  - The DM object representing the mesh

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

1019:   Level: beginner

1021: .seealso: DMPLEX, DMCreate()
1022: @*/
1023: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1024: {
1025:   PetscViewer    parentviewer = NULL;
1026:   double        *coordsIn = NULL;
1027:   GmshEntities  *entities = NULL;
1028:   GmshElement   *gmsh_elem = NULL;
1029:   PetscSection   coordSection;
1030:   Vec            coordinates;
1031:   PetscBT        periodicV = NULL, periodicC = NULL;
1032:   PetscScalar   *coords;
1033:   PetscInt       dim = 0, embedDim = -1, coordSize, c, v, d, cell, *periodicMap = NULL, *periodicMapI = NULL, *hybridMap = NULL, cMax = PETSC_DETERMINE;
1034:   PetscInt       numVertices = 0, numCells = 0, trueNumCells = 0;
1035:   int            i, shift = 1;
1036:   PetscMPIInt    rank;
1037:   PetscBool      binary, zerobase = PETSC_FALSE, usemarker = PETSC_FALSE;
1038:   PetscBool      enable_hybrid = interpolate, periodic = PETSC_TRUE;

1042:   MPI_Comm_rank(comm, &rank);
1043:   PetscObjectOptionsBegin((PetscObject)viewer);
1044:   PetscOptionsHead(PetscOptionsObject,"DMPlex Gmsh options");
1045:   PetscOptionsBool("-dm_plex_gmsh_hybrid", "Generate hybrid cell bounds", "DMPlexCreateGmsh", enable_hybrid, &enable_hybrid, NULL);
1046:   PetscOptionsBool("-dm_plex_gmsh_periodic","Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL);
1047:   PetscOptionsBool("-dm_plex_gmsh_use_marker", "Generate marker label", "DMPlexCreateGmsh", usemarker, &usemarker, NULL);
1048:   PetscOptionsBool("-dm_plex_gmsh_zero_base", "Read Gmsh file with zero base indices", "DMPlexCreateGmsh", zerobase, &zerobase, NULL);
1049:   PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", embedDim, &embedDim, NULL,PETSC_DECIDE);
1050:   PetscOptionsTail();
1051:   PetscOptionsEnd();
1052:   if (zerobase) shift = 0;

1054:   DMCreate(comm, dm);
1055:   DMSetType(*dm, DMPLEX);
1056:   PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);

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

1060:   /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1061:   if (binary) {
1062:     parentviewer = viewer;
1063:     PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1064:   }

1066:   if (!rank) {
1067:     GmshFile  gmsh_, *gmsh = &gmsh_;
1068:     char      line[PETSC_MAX_PATH_LEN];
1069:     PetscBool match;

1071:     PetscArrayzero(gmsh,1);
1072:     gmsh->viewer = viewer;
1073:     gmsh->binary = binary;

1075:     /* Read mesh format */
1076:     GmshReadSection(gmsh, line);
1077:     GmshExpect(gmsh, "$MeshFormat", line);
1078:     DMPlexCreateGmsh_ReadMeshFormat(gmsh);
1079:     GmshReadEndSection(gmsh, "$EndMeshFormat", line);

1081:     /* OPTIONAL Read physical names */
1082:     GmshReadSection(gmsh, line);
1083:     GmshMatch(gmsh,"$PhysicalNames", line, &match);
1084:     if (match) {
1085:       DMPlexCreateGmsh_ReadPhysicalNames(gmsh);
1086:       GmshReadEndSection(gmsh, "$EndPhysicalNames", line);
1087:       /* Initial read for entity section */
1088:       GmshReadSection(gmsh, line);
1089:     }

1091:     /* Read entities */
1092:     if (gmsh->fileFormat >= 40) {
1093:       GmshExpect(gmsh, "$Entities", line);
1094:       switch (gmsh->fileFormat) {
1095:       case 41: DMPlexCreateGmsh_ReadEntities_v41(gmsh, &entities); break;
1096:       default: DMPlexCreateGmsh_ReadEntities_v40(gmsh, &entities); break;
1097:       }
1098:       GmshReadEndSection(gmsh, "$EndEntities", line);
1099:       /* Initial read for nodes section */
1100:       GmshReadSection(gmsh, line);
1101:     }

1103:     /* Read nodes */
1104:     GmshExpect(gmsh, "$Nodes", line);
1105:     switch (gmsh->fileFormat) {
1106:     case 41: DMPlexCreateGmsh_ReadNodes_v41(gmsh, shift, &numVertices, &coordsIn); break;
1107:     case 40: DMPlexCreateGmsh_ReadNodes_v40(gmsh, shift, &numVertices, &coordsIn); break;
1108:     default: DMPlexCreateGmsh_ReadNodes_v22(gmsh, shift, &numVertices, &coordsIn); break;
1109:     }
1110:     GmshReadEndSection(gmsh, "$EndNodes", line);

1112:     /* Read elements */
1113:     GmshReadSection(gmsh, line);;
1114:     GmshExpect(gmsh, "$Elements", line);
1115:     switch (gmsh->fileFormat) {
1116:     case 41: DMPlexCreateGmsh_ReadElements_v41(gmsh, shift, entities, &numCells, &gmsh_elem); break;
1117:     case 40: DMPlexCreateGmsh_ReadElements_v40(gmsh, shift, entities, &numCells, &gmsh_elem); break;
1118:     default: DMPlexCreateGmsh_ReadElements_v22(gmsh, shift, &numCells, &gmsh_elem); break;
1119:     }
1120:     GmshReadEndSection(gmsh, "$EndElements", line);

1122:     /* OPTIONAL Read periodic section */
1123:     if (periodic) {
1124:       GmshReadSection(gmsh, line);
1125:       GmshMatch(gmsh, "$Periodic", line, &periodic);
1126:     }
1127:     if (periodic) {
1128:       PetscInt pVert, *periodicMapT, *aux;

1130:       PetscMalloc1(numVertices, &periodicMapT);
1131:       PetscBTCreate(numVertices, &periodicV);
1132:       for (i = 0; i < numVertices; i++) periodicMapT[i] = i;

1134:       GmshExpect(gmsh, "$Periodic", line);
1135:       switch (gmsh->fileFormat) {
1136:       case 41: DMPlexCreateGmsh_ReadPeriodic_v41(gmsh, shift, periodicMapT, periodicV); break;
1137:       default: DMPlexCreateGmsh_ReadPeriodic_v40(gmsh, shift, periodicMapT, periodicV); break;
1138:       }
1139:       GmshReadEndSection(gmsh, "$EndPeriodic", line);

1141:       /* we may have slaves of slaves */
1142:       for (i = 0; i < numVertices; i++) {
1143:         while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) {
1144:           periodicMapT[i] = periodicMapT[periodicMapT[i]];
1145:         }
1146:       }
1147:       /* periodicMap : from old to new numbering (periodic vertices excluded)
1148:          periodicMapI: from new to old numbering */
1149:       PetscMalloc1(numVertices, &periodicMap);
1150:       PetscMalloc1(numVertices, &periodicMapI);
1151:       PetscMalloc1(numVertices, &aux);
1152:       for (i = 0, pVert = 0; i < numVertices; i++) {
1153:         if (periodicMapT[i] != i) {
1154:           pVert++;
1155:         } else {
1156:           aux[i] = i - pVert;
1157:           periodicMapI[i - pVert] = i;
1158:         }
1159:       }
1160:       for (i = 0 ; i < numVertices; i++) {
1161:         periodicMap[i] = aux[periodicMapT[i]];
1162:       }
1163:       PetscFree(periodicMapT);
1164:       PetscFree(aux);
1165:       /* remove periodic vertices */
1166:       numVertices = numVertices - pVert;
1167:     }

1169:     GmshEntitiesDestroy(&entities);
1170:     PetscFree(gmsh->wbuf);
1171:     PetscFree(gmsh->sbuf);
1172:   }

1174:   if (parentviewer) {
1175:     PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1176:   }

1178:   if (!rank) {
1179:     PetscBool hybrid   = PETSC_FALSE;
1180:     PetscInt  cellType = -1;

1182:     for (trueNumCells = 0, c = 0; c < numCells; ++c) {
1183:       if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0; hybrid = PETSC_FALSE; cellType = -1;}
1184:       if (gmsh_elem[c].dim < dim) continue;
1185:       if (cellType == -1) cellType = gmsh_elem[c].cellType;
1186:       /* different cell type indicate an hybrid mesh in PLEX */
1187:       if (cellType != gmsh_elem[c].cellType) hybrid = PETSC_TRUE;
1188:       /* wedges always indicate an hybrid mesh in PLEX */
1189:       if (cellType == 6 || cellType == 13) hybrid = PETSC_TRUE;
1190:       trueNumCells++;
1191:     }
1192:     /* Renumber cells for hybrid grids */
1193:     if (hybrid && enable_hybrid) {
1194:       PetscInt hc1 = 0, hc2 = 0, *hybridCells1 = NULL, *hybridCells2 = NULL;
1195:       PetscInt cell, tn, *tp;
1196:       int      n1 = 0,n2 = 0;

1198:       PetscMalloc1(trueNumCells, &hybridCells1);
1199:       PetscMalloc1(trueNumCells, &hybridCells2);
1200:       for (cell = 0, c = 0; c < numCells; ++c) {
1201:         if (gmsh_elem[c].dim == dim) {
1202:           if (!n1) n1 = gmsh_elem[c].cellType;
1203:           else if (!n2 && n1 != gmsh_elem[c].cellType) n2 = gmsh_elem[c].cellType;

1205:           if      (gmsh_elem[c].cellType == n1) hybridCells1[hc1++] = cell;
1206:           else if (gmsh_elem[c].cellType == n2) hybridCells2[hc2++] = cell;
1207:           else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle more than 2 cell types");
1208:           cell++;
1209:         }
1210:       }

1212:       switch (n1) {
1213:       case 2: /* triangles */
1214:       case 9:
1215:         switch (n2) {
1216:         case 0: /* single type mesh */
1217:         case 3: /* quads */
1218:         case 10:
1219:           break;
1220:         default:
1221:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1222:         }
1223:         break;
1224:       case 3: /* quadrilateral */
1225:       case 10:
1226:         switch (n2) {
1227:         case 0: /* single type mesh */
1228:         case 2: /* swap since we list simplices first */
1229:         case 9:
1230:           tn  = hc1;
1231:           hc1 = hc2;
1232:           hc2 = tn;

1234:           tp           = hybridCells1;
1235:           hybridCells1 = hybridCells2;
1236:           hybridCells2 = tp;
1237:           break;
1238:         default:
1239:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1240:         }
1241:         break;
1242:       case 4: /* tetrahedra */
1243:       case 11:
1244:         switch (n2) {
1245:         case 0: /* single type mesh */
1246:         case 6: /* wedges */
1247:         case 13:
1248:           break;
1249:         default:
1250:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1251:         }
1252:         break;
1253:       case 5: /* hexahedra */
1254:       case 12:
1255:         switch (n2) {
1256:         case 0: /* single type mesh */
1257:         case 6: /* wedges */
1258:         case 13:
1259:           break;
1260:         default:
1261:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1262:         }
1263:         break;
1264:       case 6: /* wedge */
1265:       case 13:
1266:         switch (n2) {
1267:         case 0: /* single type mesh */
1268:         case 4: /* tetrahedra: swap since we list simplices first */
1269:         case 11:
1270:         case 5: /* hexahedra */
1271:         case 12:
1272:           tn  = hc1;
1273:           hc1 = hc2;
1274:           hc2 = tn;

1276:           tp           = hybridCells1;
1277:           hybridCells1 = hybridCells2;
1278:           hybridCells2 = tp;
1279:           break;
1280:         default:
1281:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1282:         }
1283:         break;
1284:       default:
1285:         SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
1286:       }
1287:       cMax = hc1;
1288:       PetscMalloc1(trueNumCells, &hybridMap);
1289:       for (cell = 0; cell < hc1; cell++) hybridMap[hybridCells1[cell]] = cell;
1290:       for (cell = 0; cell < hc2; cell++) hybridMap[hybridCells2[cell]] = cell + hc1;
1291:       PetscFree(hybridCells1);
1292:       PetscFree(hybridCells2);
1293:     }

1295:   }

1297:   /* Allocate the cell-vertex mesh */
1298:   DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
1299:   for (cell = 0, c = 0; c < numCells; ++c) {
1300:     if (gmsh_elem[c].dim == dim) {
1301:       DMPlexSetConeSize(*dm, hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].numNodes);
1302:       cell++;
1303:     }
1304:   }
1305:   DMSetUp(*dm);
1306:   /* Add cell-vertex connections */
1307:   for (cell = 0, c = 0; c < numCells; ++c) {
1308:     if (gmsh_elem[c].dim == dim) {
1309:       PetscInt pcone[8], corner;
1310:       for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1311:         const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
1312:         pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells;
1313:       }
1314:       if (dim == 3) {
1315:         /* Tetrahedra are inverted */
1316:         if (gmsh_elem[c].cellType == 4 || gmsh_elem[c].cellType == 11) {
1317:           PetscInt tmp = pcone[0];
1318:           pcone[0] = pcone[1];
1319:           pcone[1] = tmp;
1320:         }
1321:         /* Hexahedra are inverted */
1322:         if (gmsh_elem[c].cellType == 5 || gmsh_elem[c].cellType == 12) {
1323:           PetscInt tmp = pcone[1];
1324:           pcone[1] = pcone[3];
1325:           pcone[3] = tmp;
1326:         }
1327:         /* Prisms are inverted */
1328:         if (gmsh_elem[c].cellType == 6 || gmsh_elem[c].cellType == 13) {
1329:           PetscInt tmp;

1331:           tmp      = pcone[1];
1332:           pcone[1] = pcone[2];
1333:           pcone[2] = tmp;
1334:           tmp      = pcone[4];
1335:           pcone[4] = pcone[5];
1336:           pcone[5] = tmp;
1337:         }
1338:       } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */
1339:         /* quads are input to PLEX as prisms */
1340:         if (gmsh_elem[c].cellType == 3 || gmsh_elem[c].cellType == 10) {
1341:           PetscInt tmp = pcone[2];
1342:           pcone[2] = pcone[3];
1343:           pcone[3] = tmp;
1344:         }
1345:       }
1346:       DMPlexSetCone(*dm, hybridMap ? hybridMap[cell] : cell, pcone);
1347:       cell++;
1348:     }
1349:   }
1350:   MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
1351:   DMSetDimension(*dm, dim);
1352:   DMPlexSetHybridBounds(*dm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
1353:   DMPlexSymmetrize(*dm);
1354:   DMPlexStratify(*dm);
1355:   if (interpolate) {
1356:     DM idm;

1358:     DMPlexInterpolate(*dm, &idm);
1359:     DMDestroy(dm);
1360:     *dm  = idm;
1361:   }

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

1367:     DMCreateLabel(*dm, "marker");
1368:     DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
1369:     for (f = fStart; f < fEnd; ++f) {
1370:       PetscInt suppSize;

1372:       DMPlexGetSupportSize(*dm, f, &suppSize);
1373:       if (suppSize == 1) {
1374:         PetscInt *cone = NULL, coneSize, p;

1376:         DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1377:         for (p = 0; p < coneSize; p += 2) {
1378:           DMSetLabelValue(*dm, "marker", cone[p], 1);
1379:         }
1380:         DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1381:       }
1382:     }
1383:   }

1385:   if (!rank) {
1386:     PetscInt vStart, vEnd;

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

1391:       /* Create face sets */
1392:       if (interpolate && gmsh_elem[c].dim == dim-1) {
1393:         const PetscInt *join;
1394:         PetscInt        joinSize, pcone[8], corner;
1395:         /* Find the relevant facet with vertex joins */
1396:         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1397:           const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
1398:           pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart;
1399:         }
1400:         DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, pcone, &joinSize, &join);
1401:         if (joinSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for GMsh element %d (Plex cell %D)", gmsh_elem[c].id, c);
1402:         DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
1403:         DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
1404:       }

1406:       /* Create cell sets */
1407:       if (gmsh_elem[c].dim == dim) {
1408:         if (gmsh_elem[c].numTags > 0) {
1409:           DMSetLabelValue(*dm, "Cell Sets", hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].tags[0]);
1410:         }
1411:         cell++;
1412:       }

1414:       /* Create vertex sets */
1415:       if (gmsh_elem[c].dim == 0) {
1416:         if (gmsh_elem[c].numTags > 0) {
1417:           const PetscInt cc = gmsh_elem[c].nodes[0] - shift;
1418:           const PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart;
1419:           DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);
1420:         }
1421:       }
1422:     }
1423:   }

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

1471:     for (cell = 0, c = 0; c < numCells; ++c) {
1472:       PetscInt pcone[8], corner;
1473:       if (gmsh_elem[c].dim == dim) {
1474:         PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
1475:         if (PetscUnlikely(PetscBTLookup(periodicC, ucell))) {
1476:           for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1477:             pcone[corner] = gmsh_elem[c].nodes[corner] - shift;
1478:           }
1479:           if (dim == 3) {
1480:             /* Tetrahedra are inverted */
1481:             if (gmsh_elem[c].cellType == 4 || gmsh_elem[c].cellType == 11) {
1482:               PetscInt tmp = pcone[0];
1483:               pcone[0] = pcone[1];
1484:               pcone[1] = tmp;
1485:             }
1486:             /* Hexahedra are inverted */
1487:             if (gmsh_elem[c].cellType == 5 || gmsh_elem[c].cellType == 12) {
1488:               PetscInt tmp = pcone[1];
1489:               pcone[1] = pcone[3];
1490:               pcone[3] = tmp;
1491:             }
1492:             /* Prisms are inverted */
1493:             if (gmsh_elem[c].cellType == 6 || gmsh_elem[c].cellType == 13) {
1494:               PetscInt tmp;

1496:               tmp      = pcone[1];
1497:               pcone[1] = pcone[2];
1498:               pcone[2] = tmp;
1499:               tmp      = pcone[4];
1500:               pcone[4] = pcone[5];
1501:               pcone[5] = tmp;
1502:             }
1503:           } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */
1504:             /* quads are input to PLEX as prisms */
1505:             if (gmsh_elem[c].cellType == 3 || gmsh_elem[c].cellType == 10) {
1506:               PetscInt tmp = pcone[2];
1507:               pcone[2] = pcone[3];
1508:               pcone[3] = tmp;
1509:             }
1510:           }
1511:           PetscSectionGetOffset(coordSection, ucell, &off);
1512:           for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
1513:             v = pcone[corner];
1514:             for (d = 0; d < embedDim; ++d) {
1515:               coords[off++] = (PetscReal) coordsIn[v*3+d];
1516:             }
1517:           }
1518:         }
1519:         cell++;
1520:       }
1521:     }
1522:     for (v = 0; v < numVertices; ++v) {
1523:       PetscSectionGetOffset(coordSection, v + trueNumCells, &off);
1524:       for (d = 0; d < embedDim; ++d) {
1525:         coords[off+d] = (PetscReal) coordsIn[periodicMapI[v]*3+d];
1526:       }
1527:     }
1528:   } else {
1529:     for (v = 0; v < numVertices; ++v) {
1530:       for (d = 0; d < embedDim; ++d) {
1531:         coords[v*embedDim+d] = (PetscReal) coordsIn[v*3+d];
1532:       }
1533:     }
1534:   }
1535:   VecRestoreArray(coordinates, &coords);
1536:   DMSetCoordinatesLocal(*dm, coordinates);

1538:   periodic = periodicMap ? PETSC_TRUE : PETSC_FALSE;
1539:   MPI_Bcast(&periodic, 1, MPIU_BOOL, 0, comm);
1540:   DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);

1542:   PetscFree(coordsIn);
1543:   PetscFree(gmsh_elem);
1544:   PetscFree(hybridMap);
1545:   PetscFree(periodicMap);
1546:   PetscFree(periodicMapI);
1547:   PetscBTDestroy(&periodicV);
1548:   PetscBTDestroy(&periodicC);
1549:   VecDestroy(&coordinates);

1551:   PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
1552:   return(0);
1553: }