Actual source code: plexgmsh.c

petsc-3.9.2 2018-05-20
Report Typos and Errors
  1: #define PETSCDM_DLL
  2:  #include <petsc/private/dmpleximpl.h>

  4: /*@C
  5:   DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file

  7: + comm        - The MPI communicator
  8: . filename    - Name of the Gmsh file
  9: - interpolate - Create faces and edges in the mesh

 11:   Output Parameter:
 12: . dm  - The DM object representing the mesh

 14:   Level: beginner

 16: .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
 17: @*/
 18: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
 19: {
 20:   PetscViewer     viewer;
 21:   PetscMPIInt     rank;
 22:   int             fileType;
 23:   PetscViewerType vtype;
 24:   PetscErrorCode  ierr;

 27:   MPI_Comm_rank(comm, &rank);

 29:   /* Determine Gmsh file type (ASCII or binary) from file header */
 30:   if (!rank) {
 31:     PetscViewer vheader;
 32:     char        line[PETSC_MAX_PATH_LEN];
 33:     PetscBool   match;
 34:     int         snum;
 35:     float       version;

 37:     PetscViewerCreate(PETSC_COMM_SELF, &vheader);
 38:     PetscViewerSetType(vheader, PETSCVIEWERASCII);
 39:     PetscViewerFileSetMode(vheader, FILE_MODE_READ);
 40:     PetscViewerFileSetName(vheader, filename);
 41:     /* Read only the first two lines of the Gmsh file */
 42:     PetscViewerRead(vheader, line, 1, NULL, PETSC_STRING);
 43:     PetscStrncmp(line, "$MeshFormat", sizeof(line), &match);
 44:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
 45:     PetscViewerRead(vheader, line, 2, NULL, PETSC_STRING);
 46:     snum = sscanf(line, "%f %d", &version, &fileType);
 47:     if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
 48:     if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
 49:     PetscViewerDestroy(&vheader);
 50:   }
 51:   MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);
 52:   vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;

 54:   /* Create appropriate viewer and build plex */
 55:   PetscViewerCreate(comm, &viewer);
 56:   PetscViewerSetType(viewer, vtype);
 57:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);
 58:   PetscViewerFileSetName(viewer, filename);
 59:   DMPlexCreateGmsh(comm, viewer, interpolate, dm);
 60:   PetscViewerDestroy(&viewer);
 61:   return(0);
 62: }

 64: static PetscErrorCode DMPlexCreateGmsh_ReadNodes(PetscViewer viewer, PetscBool binary, PetscBool byteSwap, int shift, int numVertices, double **coordinates)
 65: {
 66:   int            v,nid;

 70:   PetscMalloc1(numVertices*3, coordinates);
 71:   for (v = 0; v < numVertices; ++v) {
 72:     double *xyz = *coordinates + v*3;
 73:     PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
 74:     if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
 75:     if (nid != v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %d", nid, v+shift);
 76:     PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
 77:     if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
 78:   }
 79:   return(0);
 80: }

 82: static PetscErrorCode DMPlexCreateGmsh_ReadElements(PetscViewer viewer, PetscBool binary, PetscBool byteSwap, PETSC_UNUSED int shift, int numCells, GmshElement **gmsh_elems)
 83: {
 84:   GmshElement   *elements;
 85:   int            i, c, p, ibuf[1+4+512];
 86:   int            cellType, dim, numNodes, numNodesIgnore, numElem, numTags;

 90:   PetscMalloc1(numCells, &elements);
 91:   for (c = 0; c < numCells;) {
 92:     PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
 93:     if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 3);}
 94:     if (binary) {
 95:       cellType = ibuf[0];
 96:       numElem = ibuf[1];
 97:       numTags = ibuf[2];
 98:     } else {
 99:       elements[c].id = ibuf[0];
100:       cellType = ibuf[1];
101:       numTags = ibuf[2];
102:       numElem = 1;
103:     }
104:     /* http://gmsh.info/doc/texinfo/gmsh.html#MSH-ASCII-file-format */
105:     numNodesIgnore = 0;
106:     switch (cellType) {
107:     case 1: /* 2-node line */
108:       dim = 1;
109:       numNodes = 2;
110:       break;
111:     case 2: /* 3-node triangle */
112:       dim = 2;
113:       numNodes = 3;
114:       break;
115:     case 3: /* 4-node quadrangle */
116:       dim = 2;
117:       numNodes = 4;
118:       break;
119:     case 4: /* 4-node tetrahedron */
120:       dim  = 3;
121:       numNodes = 4;
122:       break;
123:     case 5: /* 8-node hexahedron */
124:       dim = 3;
125:       numNodes = 8;
126:       break;
127:     case 8: /* 3-node 2nd order line */
128:       dim = 1;
129:       numNodes = 2;
130:       numNodesIgnore = 1;
131:       break;
132:     case 9: /* 6-node 2nd order triangle */
133:       dim = 2;
134:       numNodes = 3;
135:       numNodesIgnore = 3;
136:       break;
137:     case 15: /* 1-node vertex */
138:       dim = 0;
139:       numNodes = 1;
140:       break;
141:     case 6: /* 6-node prism */
142:     case 7: /* 5-node pyramid */
143:     case 10: /* 9-node 2nd order quadrangle */
144:     case 11: /* 10-node 2nd order tetrahedron */
145:     case 12: /* 27-node 2nd order hexhedron */
146:     case 13: /* 19-node 2nd order prism */
147:     case 14: /* 14-node 2nd order pyramid */
148:     default:
149:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
150:     }
151:     if (binary) {
152:       const int nint = 1 + numTags + numNodes + numNodesIgnore;
153:       /* Loop over element blocks */
154:       for (i = 0; i < numElem; ++i, ++c) {
155:         PetscViewerRead(viewer, ibuf, nint, NULL, PETSC_ENUM);
156:         if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, nint);}
157:         elements[c].dim = dim;
158:         elements[c].numNodes = numNodes;
159:         elements[c].numTags = numTags;
160:         elements[c].id = ibuf[0];
161:         for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p];
162:         for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
163:       }
164:     } else {
165:       elements[c].dim = dim;
166:       elements[c].numNodes = numNodes;
167:       elements[c].numTags = numTags;
168:       PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);
169:       PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);
170:       PetscViewerRead(viewer, ibuf, numNodesIgnore, NULL, PETSC_ENUM);
171:       c++;
172:     }
173:   }
174:   *gmsh_elems = elements;
175:   return(0);
176: }

178: /*@
179:   DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer

181:   Collective on comm

183:   Input Parameters:
184: + comm  - The MPI communicator
185: . viewer - The Viewer associated with a Gmsh file
186: - interpolate - Create faces and edges in the mesh

188:   Output Parameter:
189: . dm  - The DM object representing the mesh

191:   Note: http://www.geuz.org/gmsh/doc/texinfo/#MSH-ASCII-file-format
192:   and http://www.geuz.org/gmsh/doc/texinfo/#MSH-binary-file-format

194:   Level: beginner

196: .keywords: mesh,Gmsh
197: .seealso: DMPLEX, DMCreate()
198: @*/
199: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
200: {
201:   PetscViewer    parentviewer = NULL;
202:   double        *coordsIn = NULL;
203:   GmshElement   *gmsh_elem = NULL;
204:   PetscSection   coordSection;
205:   Vec            coordinates;
206:   PetscBT        periodicV = NULL, periodicC = NULL;
207:   PetscScalar   *coords;
208:   PetscInt       dim = 0, embedDim, coordSize, c, v, d, r, cell, *periodicMap = NULL, *periodicMapI = NULL;
209:   int            i, numVertices = 0, numCells = 0, trueNumCells = 0, numRegions = 0, snum, shift = 1;
210:   PetscMPIInt    rank;
211:   char           line[PETSC_MAX_PATH_LEN];
212:   PetscBool      binary, byteSwap = PETSC_FALSE, zerobase = PETSC_FALSE, isbd = PETSC_FALSE, periodic = PETSC_FALSE, usemarker = PETSC_FALSE;

216:   MPI_Comm_rank(comm, &rank);
217:   PetscOptionsGetBool(NULL, NULL, "-dm_plex_gmsh_periodic", &periodic, NULL);
218:   PetscOptionsGetBool(NULL, NULL, "-dm_plex_gmsh_use_marker", &usemarker, NULL);
219:   PetscOptionsGetBool(NULL, NULL, "-gmsh_zero_base", &zerobase, NULL);
220:   PetscOptionsGetBool(NULL, NULL, "-gmsh_bd", &isbd, NULL);
221:   if (zerobase) shift = 0;

223:   DMCreate(comm, dm);
224:   DMSetType(*dm, DMPLEX);
225:   PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);

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

229:   /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
230:   if (binary) {
231:     parentviewer = viewer;
232:     PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
233:   }

235:   if (!rank) {
236:     PetscBool match;
237:     int       fileType, dataSize;
238:     float     version;

240:     /* Read header */
241:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
242:     PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);
243:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
244:     PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
245:     snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
246:     if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
247:     if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
248:     if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
249:     if (binary) {
250:       int checkInt;
251:       PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);
252:       if (checkInt != 1) {
253:         PetscByteSwap(&checkInt, PETSC_ENUM, 1);
254:         if (checkInt == 1) byteSwap = PETSC_TRUE;
255:         else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType);
256:       }
257:     } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType);
258:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
259:     PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);
260:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");

262:     /* OPTIONAL Read physical names */
263:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
264:     PetscStrncmp(line, "$PhysicalNames", PETSC_MAX_PATH_LEN, &match);
265:     if (match) {
266:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
267:       snum = sscanf(line, "%d", &numRegions);
268:       if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
269:       for (r = 0; r < numRegions; ++r) {
270:         PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
271:       }
272:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
273:       PetscStrncmp(line, "$EndPhysicalNames", PETSC_MAX_PATH_LEN, &match);
274:       if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
275:       /* Initial read for vertex section */
276:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
277:     }

279:     /* Read vertices */
280:     PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);
281:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
282:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
283:     snum = sscanf(line, "%d", &numVertices);
284:     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
285:     DMPlexCreateGmsh_ReadNodes(viewer, binary, byteSwap, shift, numVertices, &coordsIn);
286:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
287:     PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);
288:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");

290:     /* Read cells */
291:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
292:     PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);
293:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
294:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
295:     snum = sscanf(line, "%d", &numCells);
296:     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
297:     /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
298:        file contents multiple times to figure out the true number of cells and facets
299:        in the given mesh. To make this more efficient we read the file contents only
300:        once and store them in memory, while determining the true number of cells. */
301:     DMPlexCreateGmsh_ReadElements(viewer, binary, byteSwap, shift, numCells, &gmsh_elem);
302:     for (trueNumCells=0, c = 0; c < numCells; ++c) {
303:       if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;}
304:       if (gmsh_elem[c].dim == dim) trueNumCells++;
305:     }
306:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
307:     PetscStrncmp(line, "$EndElements", 12, &match);
308:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");

310:     /* OPTIONAL Read periodic section */
311:     if (periodic) {
312:       PetscInt pVert, *periodicMapT, *aux;
313:       int      numPeriodic;

315:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
316:       PetscStrncmp(line, "$Periodic", 9, &match);
317:       if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
318:       PetscMalloc1(numVertices, &periodicMapT);
319:       PetscBTCreate(numVertices, &periodicV);
320:       for (i = 0; i < numVertices; i++) periodicMapT[i] = i;
321:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
322:       snum = sscanf(line, "%d", &numPeriodic);
323:       if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
324:       for (i = 0; i < numPeriodic; i++) {
325:         int j, edim = -1, slaveTag = -1, masterTag = -1, nNodes, slaveNode, masterNode;

327:         PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
328:         snum = sscanf(line, "%d %d %d", &edim, &slaveTag, &masterTag); /* slaveTag and masterTag are unused */
329:         if (snum != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
330:         PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
331:         snum = sscanf(line, "%d", &nNodes);
332:         if (snum != 1) { /* discard tranformation and try again */
333:           PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
334:           PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
335:           snum = sscanf(line, "%d", &nNodes);
336:           if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
337:         }
338:         for (j = 0; j < nNodes; j++) {
339:           PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
340:           snum = sscanf(line, "%d %d", &slaveNode, &masterNode);
341:           if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
342:           periodicMapT[slaveNode - shift] = masterNode - shift;
343:           PetscBTSet(periodicV, slaveNode - shift);
344:           PetscBTSet(periodicV, masterNode - shift);
345:         }
346:       }
347:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
348:       PetscStrncmp(line, "$EndPeriodic", 12, &match);
349:       if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
350:       /* we may have slaves of slaves */
351:       for (i = 0; i < numVertices; i++) {
352:         while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) {
353:           periodicMapT[i] = periodicMapT[periodicMapT[i]];
354:         }
355:       }
356:       /* periodicMap : from old to new numbering (periodic vertices excluded)
357:          periodicMapI: from new to old numbering */
358:       PetscMalloc1(numVertices, &periodicMap);
359:       PetscMalloc1(numVertices, &periodicMapI);
360:       PetscMalloc1(numVertices, &aux);
361:       for (i = 0, pVert = 0; i < numVertices; i++) {
362:         if (periodicMapT[i] != i) {
363:           pVert++;
364:         } else {
365:           aux[i] = i - pVert;
366:           periodicMapI[i - pVert] = i;
367:         }
368:       }
369:       for (i = 0 ; i < numVertices; i++) {
370:         periodicMap[i] = aux[periodicMapT[i]];
371:       }
372:       PetscFree(periodicMapT);
373:       PetscFree(aux);
374:       /* remove periodic vertices */
375:       numVertices = numVertices - pVert;
376:     }
377:   }

379:   if (parentviewer) {
380:     PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
381:   }

383:   /* Allocate the cell-vertex mesh */
384:   DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
385:   for (cell = 0, c = 0; c < numCells; ++c) {
386:     if (gmsh_elem[c].dim == dim) {
387:       DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);
388:       cell++;
389:     }
390:   }
391:   DMSetUp(*dm);
392:   /* Add cell-vertex connections */
393:   for (cell = 0, c = 0; c < numCells; ++c) {
394:     if (gmsh_elem[c].dim == dim) {
395:       PetscInt pcone[8], corner;
396:       for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
397:         const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
398:         pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells;
399:       }
400:       if (dim == 3) {
401:         /* Tetrahedra are inverted */
402:         if (gmsh_elem[c].numNodes == 4) {
403:           PetscInt tmp = pcone[0];
404:           pcone[0] = pcone[1];
405:           pcone[1] = tmp;
406:         }
407:         /* Hexahedra are inverted */
408:         if (gmsh_elem[c].numNodes == 8) {
409:           PetscInt tmp = pcone[1];
410:           pcone[1] = pcone[3];
411:           pcone[3] = tmp;
412:         }
413:       }
414:       DMPlexSetCone(*dm, cell, pcone);
415:       cell++;
416:     }
417:   }
418:   MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
419:   DMSetDimension(*dm, dim);
420:   DMPlexSymmetrize(*dm);
421:   DMPlexStratify(*dm);
422:   if (interpolate) {
423:     DM idm;

425:     DMPlexInterpolate(*dm, &idm);
426:     DMDestroy(dm);
427:     *dm  = idm;
428:   }

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

434:     DMCreateLabel(*dm, "marker");
435:     DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
436:     for (f = fStart; f < fEnd; ++f) {
437:       PetscInt suppSize;

439:       DMPlexGetSupportSize(*dm, f, &suppSize);
440:       if (suppSize == 1) {
441:         PetscInt *cone = NULL, coneSize, p;

443:         DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
444:         for (p = 0; p < coneSize; p += 2) {
445:           DMSetLabelValue(*dm, "marker", cone[p], 1);
446:         }
447:         DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
448:       }
449:     }
450:   }

452:   if (!rank) {
453:     PetscInt vStart, vEnd;

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

458:       /* Create face sets */
459:       if (interpolate && gmsh_elem[c].dim == dim-1) {
460:         const PetscInt *join;
461:         PetscInt        joinSize, pcone[8], corner;
462:         /* Find the relevant facet with vertex joins */
463:         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
464:           const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
465:           pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart;
466:         }
467:         DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, pcone, &joinSize, &join);
468:         if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id);
469:         DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
470:         DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
471:       }

473:       /* Create cell sets */
474:       if (gmsh_elem[c].dim == dim) {
475:         if (gmsh_elem[c].numTags > 0) {
476:           DMSetLabelValue(*dm, "Cell Sets", cell, gmsh_elem[c].tags[0]);
477:         }
478:         cell++;
479:       }

481:       /* Create vertex sets */
482:       if (gmsh_elem[c].dim == 0) {
483:         if (gmsh_elem[c].numTags > 0) {
484:           const PetscInt cc = gmsh_elem[c].nodes[0] - shift;
485:           const PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart;
486:           DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);
487:         }
488:       }
489:     }
490:   }

492:   /* Create coordinates */
493:   embedDim = isbd ? dim+1 : dim;
494:   DMGetCoordinateSection(*dm, &coordSection);
495:   PetscSectionSetNumFields(coordSection, 1);
496:   PetscSectionSetFieldComponents(coordSection, 0, embedDim);
497:   if (periodicMap) { /* we need to localize coordinates on cells */
498:     PetscSectionSetChart(coordSection, 0, trueNumCells + numVertices);
499:   } else {
500:     PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);
501:   }
502:   for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
503:     PetscSectionSetDof(coordSection, v, embedDim);
504:     PetscSectionSetFieldDof(coordSection, v, 0, embedDim);
505:   }
506:   if (periodicMap) {
507:     PetscBTCreate(trueNumCells, &periodicC);
508:     for (cell = 0, c = 0; c < numCells; ++c) {
509:       if (gmsh_elem[c].dim == dim) {
510:         PetscInt  corner;
511:         PetscBool pc = PETSC_FALSE;
512:         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
513:           pc = (PetscBool)(pc || PetscBTLookup(periodicV, gmsh_elem[c].nodes[corner] - shift));
514:         }
515:         if (pc) {
516:           PetscInt dof = gmsh_elem[c].numNodes*embedDim;
517:           PetscBTSet(periodicC, cell);
518:           PetscSectionSetDof(coordSection, cell, dof);
519:           PetscSectionSetFieldDof(coordSection, cell, 0, dof);
520:         }
521:         cell++;
522:       }
523:     }
524:   }
525:   PetscSectionSetUp(coordSection);
526:   PetscSectionGetStorageSize(coordSection, &coordSize);
527:   VecCreate(PETSC_COMM_SELF, &coordinates);
528:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
529:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
530:   VecSetBlockSize(coordinates, embedDim);
531:   VecSetType(coordinates, VECSTANDARD);
532:   VecGetArray(coordinates, &coords);
533:   if (periodicMap) {
534:     PetscInt off;

536:     for (cell = 0, c = 0; c < numCells; ++c) {
537:       PetscInt pcone[8], corner;
538:       if (gmsh_elem[c].dim == dim) {
539:         if (PetscUnlikely(PetscBTLookup(periodicC, cell))) {
540:           for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
541:             pcone[corner] = gmsh_elem[c].nodes[corner] - shift;
542:           }
543:           if (dim == 3) {
544:             /* Tetrahedra are inverted */
545:             if (gmsh_elem[c].numNodes == 4) {
546:               PetscInt tmp = pcone[0];
547:               pcone[0] = pcone[1];
548:               pcone[1] = tmp;
549:             }
550:             /* Hexahedra are inverted */
551:             if (gmsh_elem[c].numNodes == 8) {
552:               PetscInt tmp = pcone[1];
553:               pcone[1] = pcone[3];
554:               pcone[3] = tmp;
555:             }
556:           }
557:           PetscSectionGetOffset(coordSection, cell, &off);
558:           for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
559:             v = pcone[corner];
560:             for (d = 0; d < embedDim; ++d) {
561:               coords[off++] = (PetscReal) coordsIn[v*3+d];
562:             }
563:           }
564:         }
565:         cell++;
566:       }
567:     }
568:     for (v = 0; v < numVertices; ++v) {
569:       PetscSectionGetOffset(coordSection, v + trueNumCells, &off);
570:       for (d = 0; d < embedDim; ++d) {
571:         coords[off+d] = (PetscReal) coordsIn[periodicMapI[v]*3+d];
572:       }
573:     }
574:   } else {
575:     for (v = 0; v < numVertices; ++v) {
576:       for (d = 0; d < embedDim; ++d) {
577:         coords[v*embedDim+d] = (PetscReal) coordsIn[v*3+d];
578:       }
579:     }
580:   }
581:   VecRestoreArray(coordinates, &coords);
582:   DMSetCoordinatesLocal(*dm, coordinates);
583:   DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);

585:   PetscFree(coordsIn);
586:   PetscFree(gmsh_elem);
587:   PetscFree(periodicMap);
588:   PetscFree(periodicMapI);
589:   PetscBTDestroy(&periodicV);
590:   PetscBTDestroy(&periodicC);
591:   VecDestroy(&coordinates);

593:   PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
594:   return(0);
595: }