Actual source code: plexgmsh.c

petsc-3.8.3 2017-12-09
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, vheader;
 21:   PetscMPIInt     rank;
 22:   PetscViewerType vtype;
 23:   char            line[PETSC_MAX_PATH_LEN];
 24:   int             snum;
 25:   PetscBool       match;
 26:   int             fT;
 27:   PetscInt        fileType;
 28:   float           version;
 29:   PetscErrorCode  ierr;

 32:   MPI_Comm_rank(comm, &rank);
 33:   /* Determine Gmsh file type (ASCII or binary) from file header */
 34:   PetscViewerCreate(comm, &vheader);
 35:   PetscViewerSetType(vheader, PETSCVIEWERASCII);
 36:   PetscViewerFileSetMode(vheader, FILE_MODE_READ);
 37:   PetscViewerFileSetName(vheader, filename);
 38:   if (!rank) {
 39:     /* Read only the first two lines of the Gmsh file */
 40:     PetscViewerRead(vheader, line, 1, NULL, PETSC_STRING);
 41:     PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);
 42:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
 43:     PetscViewerRead(vheader, line, 2, NULL, PETSC_STRING);
 44:     snum = sscanf(line, "%f %d", &version, &fT);
 45:     fileType = (PetscInt) fT;
 46:     if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
 47:     if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
 48:   }
 49:   MPI_Bcast(&fileType, 1, MPIU_INT, 0, comm);
 50:   /* Create appropriate viewer and build plex */
 51:   if (fileType == 0) vtype = PETSCVIEWERASCII;
 52:   else vtype = PETSCVIEWERBINARY;
 53:   PetscViewerCreate(comm, &viewer);
 54:   PetscViewerSetType(viewer, vtype);
 55:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);
 56:   PetscViewerFileSetName(viewer, filename);
 57:   DMPlexCreateGmsh(comm, viewer, interpolate, dm);
 58:   PetscViewerDestroy(&viewer);
 59:   PetscViewerDestroy(&vheader);
 60:   return(0);
 61: }

 63: static PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscInt numCells, PetscBool binary, PetscBool byteSwap, GmshElement **gmsh_elems)
 64: {
 65:   PetscInt       c, p;
 66:   GmshElement   *elements;
 67:   int            i, cellType, dim, numNodes, numNodesIgnore, numElem, numTags;
 68:   PetscInt       pibuf[64];
 69:   int            ibuf[16];

 73:   PetscMalloc1(numCells, &elements);
 74:   for (c = 0; c < numCells;) {
 75:     PetscViewerRead(viewer, &ibuf, 3, NULL, PETSC_ENUM);
 76:     if (byteSwap) PetscByteSwap(&ibuf, PETSC_ENUM, 3);
 77:     if (binary) {
 78:       cellType = ibuf[0];
 79:       numElem = ibuf[1];
 80:       numTags = ibuf[2];
 81:     } else {
 82:       elements[c].id = ibuf[0];
 83:       cellType = ibuf[1];
 84:       numTags = ibuf[2];
 85:       numElem = 1;
 86:     }
 87:     /* http://gmsh.info/doc/texinfo/gmsh.html#MSH-ASCII-file-format */
 88:     numNodesIgnore = 0;
 89:     switch (cellType) {
 90:     case 1: /* 2-node line */
 91:       dim = 1;
 92:       numNodes = 2;
 93:       break;
 94:     case 2: /* 3-node triangle */
 95:       dim = 2;
 96:       numNodes = 3;
 97:       break;
 98:     case 3: /* 4-node quadrangle */
 99:       dim = 2;
100:       numNodes = 4;
101:       break;
102:     case 4: /* 4-node tetrahedron */
103:       dim  = 3;
104:       numNodes = 4;
105:       break;
106:     case 5: /* 8-node hexahedron */
107:       dim = 3;
108:       numNodes = 8;
109:       break;
110:     case 8: /* 3-node 2nd order line */
111:       dim = 1;
112:       numNodes = 2;
113:       numNodesIgnore = 1;
114:       break;
115:     case 9: /* 6-node 2nd order triangle */
116:       dim = 2;
117:       numNodes = 3;
118:       numNodesIgnore = 3;
119:       break;
120:     case 15: /* 1-node vertex */
121:       dim = 0;
122:       numNodes = 1;
123:       break;
124:     case 6: /* 6-node prism */
125:     case 7: /* 5-node pyramid */
126:     case 10: /* 9-node 2nd order quadrangle */
127:     case 11: /* 10-node 2nd order tetrahedron */
128:     case 12: /* 27-node 2nd order hexhedron */
129:     case 13: /* 19-node 2nd order prism */
130:     case 14: /* 14-node 2nd order pyramid */
131:     default:
132:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
133:     }
134:     if (binary) {
135:       const PetscInt nint = numNodes + numTags + 1 + numNodesIgnore;
136:       for (i = 0; i < numElem; ++i, ++c) {
137:         /* Loop over inner binary element block */
138:         elements[c].dim = dim;
139:         elements[c].numNodes = numNodes;
140:         elements[c].numTags = numTags;

142:         PetscViewerRead(viewer, &ibuf, nint, NULL, PETSC_ENUM);
143:         if (byteSwap) PetscByteSwap( &ibuf, PETSC_ENUM, nint);
144:         elements[c].id = ibuf[0];
145:         for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p];
146:         for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
147:       }
148:     } else {
149:       elements[c].dim = dim;
150:       elements[c].numNodes = numNodes;
151:       elements[c].numTags = numTags;
152:       PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);
153:       PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);
154:       PetscViewerRead(viewer, pibuf, numNodesIgnore, NULL, PETSC_ENUM);
155:       c++;
156:     }
157:   }
158:   *gmsh_elems = elements;
159:   return(0);
160: }

162: /*@
163:   DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer

165:   Collective on comm

167:   Input Parameters:
168: + comm  - The MPI communicator
169: . viewer - The Viewer associated with a Gmsh file
170: - interpolate - Create faces and edges in the mesh

172:   Output Parameter:
173: . dm  - The DM object representing the mesh

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

178:   Level: beginner

180: .keywords: mesh,Gmsh
181: .seealso: DMPLEX, DMCreate()
182: @*/
183: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
184: {
185:   PetscViewerType vtype;
186:   GmshElement   *gmsh_elem = NULL;
187:   PetscSection   coordSection;
188:   Vec            coordinates;
189:   PetscBT        periodicV = NULL, periodicC = NULL;
190:   PetscScalar   *coords;
191:   PetscReal     *coordsIn = NULL;
192:   PetscInt       dim = 0, embedDim, coordSize, c, v, d, r, cell, *periodicMap = NULL, *periodicMapI = NULL;
193:   int            i, numVertices = 0, numCells = 0, trueNumCells = 0, numRegions = 0, snum, shift = 1;
194:   PetscMPIInt    num_proc, rank;
195:   char           line[PETSC_MAX_PATH_LEN];
196:   PetscBool      zerobase = PETSC_FALSE, isbd = PETSC_FALSE, match, binary, bswap = PETSC_FALSE, periodic = PETSC_FALSE, usemarker = PETSC_FALSE;

200:   PetscOptionsGetBool(NULL, NULL, "-gmsh_zero_base", &zerobase, NULL);
201:   if (zerobase) shift = 0;
202:   MPI_Comm_rank(comm, &rank);
203:   MPI_Comm_size(comm, &num_proc);
204:   DMCreate(comm, dm);
205:   DMSetType(*dm, DMPLEX);
206:   PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);
207:   PetscViewerGetType(viewer, &vtype);
208:   PetscStrcmp(vtype, PETSCVIEWERBINARY, &binary);
209:   PetscOptionsGetBool(NULL,NULL,"-dm_plex_gmsh_periodic",&periodic,NULL);
210:   PetscOptionsGetBool(NULL,NULL,"-dm_plex_gmsh_usemarker",&usemarker,NULL);
211:   if (!rank || binary) {
212:     PetscBool match;
213:     int       fileType, dataSize;
214:     float     version;

216:     /* Read header */
217:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
218:     PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);
219:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
220:     PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
221:     snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
222:     if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
223:     if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
224:     if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
225:     if (binary) {
226:       int checkInt;
227:       PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);
228:       if (checkInt != 1) {
229:         PetscByteSwap(&checkInt, PETSC_ENUM, 1);
230:         if (checkInt == 1) bswap = PETSC_TRUE;
231:         else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType);
232:       }
233:     } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType);
234:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
235:     PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);
236:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
237:     /* OPTIONAL Read physical names */
238:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
239:     PetscStrncmp(line, "$PhysicalNames", PETSC_MAX_PATH_LEN, &match);
240:     if (match) {
241:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
242:       snum = sscanf(line, "%d", &numRegions);
243:       if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
244:       for (r = 0; r < numRegions; ++r) {PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);}
245:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
246:       PetscStrncmp(line, "$EndPhysicalNames", PETSC_MAX_PATH_LEN, &match);
247:       if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
248:       /* Initial read for vertex section */
249:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
250:     }
251:     /* Read vertices */
252:     PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);
253:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
254:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
255:     snum = sscanf(line, "%d", &numVertices);
256:     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
257:     PetscMalloc1(numVertices*3, &coordsIn);
258:     if (binary) {
259:       size_t   doubleSize, intSize;
260:       PetscInt elementSize;
261:       char     *buffer;
262:       double   *baseptr;
263:       PetscDataTypeGetSize(PETSC_ENUM, &intSize);
264:       PetscDataTypeGetSize(PETSC_DOUBLE, &doubleSize);
265:       elementSize = (intSize + 3*doubleSize);
266:       PetscMalloc1(elementSize*numVertices, &buffer);
267:       PetscViewerRead(viewer, buffer, elementSize*numVertices, NULL, PETSC_CHAR);
268:       if (bswap) PetscByteSwap(buffer, PETSC_CHAR, elementSize*numVertices);
269:       for (v = 0; v < numVertices; ++v) {
270:         baseptr = ((double*)(buffer+v*elementSize+intSize));
271:         coordsIn[v*3+0] = (PetscReal) baseptr[0];
272:         coordsIn[v*3+1] = (PetscReal) baseptr[1];
273:         coordsIn[v*3+2] = (PetscReal) baseptr[2];
274:       }
275:       PetscFree(buffer);
276:     } else {
277:       for (v = 0; v < numVertices; ++v) {
278:         PetscViewerRead(viewer, &i, 1, NULL, PETSC_ENUM);
279:         PetscViewerRead(viewer, &(coordsIn[v*3]), 3, NULL, PETSC_DOUBLE);
280:         if (i != (int) v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+shift);
281:       }
282:     }
283:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
284:     PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);
285:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
286:     /* Read cells */
287:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
288:     PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);
289:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
290:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
291:     snum = sscanf(line, "%d", &numCells);
292:     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
293:   }

295:   if (!rank || binary) {
296:     /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
297:        file contents multiple times to figure out the true number of cells and facets
298:        in the given mesh. To make this more efficient we read the file contents only
299:        once and store them in memory, while determining the true number of cells. */
300:     DMPlexCreateGmsh_ReadElement(viewer, numCells, binary, bswap, &gmsh_elem);
301:     for (trueNumCells=0, c = 0; c < numCells; ++c) {
302:       if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;}
303:       if (gmsh_elem[c].dim == dim) trueNumCells++;
304:     }
305:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
306:     PetscStrncmp(line, "$EndElements", 12, &match);
307:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
308:     if (periodic) {
309:       PetscInt pVert, *periodicMapT, *aux;
310:       int      numPeriodic;

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

324:         PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
325:         snum = sscanf(line, "%d %d %d", &edim, &slaveTag, &masterTag); /* slaveTag and masterTag are unused */
326:         if (snum != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
327:         /* Type of tranformation, discarded */
328:         PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
329:         PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
330:         snum = sscanf(line, "%d", &nNodes);
331:         if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
332:         for (j = 0; j < nNodes; j++) {
333:           PetscInt slaveNode, masterNode;

335:           PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
336:           snum = sscanf(line, "%" PetscInt_FMT " %" PetscInt_FMT, &slaveNode, &masterNode);
337:           if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
338:           periodicMapT[slaveNode - shift] = masterNode - shift;
339:           PetscBTSet(periodicV, slaveNode - shift);
340:           PetscBTSet(periodicV, masterNode - shift);
341:         }
342:       }
343:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
344:       PetscStrncmp(line, "$EndPeriodic", 12, &match);
345:       if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
346:       /* we may have slaves of slaves */
347:       for (i = 0; i < numVertices; i++) {
348:         while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) {
349:           periodicMapT[i] = periodicMapT[periodicMapT[i]];
350:         }
351:       }
352:       /* periodicMap : from old to new numbering (periodic vertices excluded)
353:          periodicMapI: from new to old numbering */
354:       PetscMalloc1(numVertices, &periodicMap);
355:       PetscMalloc1(numVertices, &periodicMapI);
356:       PetscMalloc1(numVertices, &aux);
357:       for (i = 0, pVert = 0; i < numVertices; i++) {
358:         if (periodicMapT[i] != i) {
359:           pVert++;
360:         } else {
361:           aux[i] = i - pVert;
362:           periodicMapI[i - pVert] = i;
363:         }
364:       }
365:       for (i = 0 ; i < numVertices; i++) {
366:         periodicMap[i] = aux[periodicMapT[i]];
367:       }
368:       PetscFree(periodicMapT);
369:       PetscFree(aux);
370:       /* remove periodic vertices */
371:       numVertices = numVertices - pVert;
372:     }
373:   }
374:   /* For binary we read on all ranks, but only build the plex on rank 0 */
375:   if (binary && rank) {trueNumCells = 0; numVertices = 0;};
376:   /* Allocate the cell-vertex mesh */
377:   DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
378:   if (!rank) {
379:     for (cell = 0, c = 0; c < numCells; ++c) {
380:       if (gmsh_elem[c].dim == dim) {
381:         DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);
382:         cell++;
383:       }
384:     }
385:   }
386:   DMSetUp(*dm);
387:   /* Add cell-vertex connections */
388:   if (!rank) {
389:     PetscInt pcone[8], corner;

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

423:     DMPlexInterpolate(*dm, &idm);
424:     DMDestroy(dm);
425:     *dm  = idm;
426:   }

428:   if (!rank && usemarker) {
429:     PetscInt f, fStart, fEnd;

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

436:       DMPlexGetSupportSize(*dm, f, &suppSize);
437:       if (suppSize == 1) {
438:         PetscInt *cone = NULL, coneSize, p;

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

449:   if (!rank) {
450:     /* Apply boundary IDs by finding the relevant facets with vertex joins */
451:     PetscInt pcone[8], corner, vStart, vEnd;

453:     DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
454:     for (c = 0; c < numCells; ++c) {
455:       if (gmsh_elem[c].dim == dim-1) {
456:         const PetscInt *join;
457:         PetscInt        joinSize;

459:         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
460:           const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
461:           pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart;
462:         }
463:         DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
464:         if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id);
465:         DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
466:         DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
467:       }
468:     }

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

480:     /* Create vertex sets */
481:     for (c = 0; c < numCells; ++c) {
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:           PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart;
486:           DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);
487:         }
488:       }
489:     }
490:   }

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

536:     if (periodicMap) {
537:       PetscInt off;

539:       for (cell = 0, c = 0; c < numCells; ++c) {
540:         PetscInt pcone[8], corner;
541:         if (gmsh_elem[c].dim == dim) {
542:           if (PetscUnlikely(PetscBTLookup(periodicC,cell))) {
543:             for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
544:               pcone[corner] = gmsh_elem[c].nodes[corner] - shift;
545:             }
546:             if (dim == 3) {
547:               /* Tetrahedra are inverted */
548:               if (gmsh_elem[c].numNodes == 4) {
549:                 PetscInt tmp = pcone[0];
550:                 pcone[0] = pcone[1];
551:                 pcone[1] = tmp;
552:               }
553:               /* Hexahedra are inverted */
554:               if (gmsh_elem[c].numNodes == 8) {
555:                 PetscInt tmp = pcone[1];
556:                 pcone[1] = pcone[3];
557:                 pcone[3] = tmp;
558:               }
559:             }
560:             PetscSectionGetOffset(coordSection, cell, &off);
561:             for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
562:               PetscInt v = pcone[corner];
563:               for (d = 0; d < embedDim; ++d) {
564:                 coords[off++] = coordsIn[v*3+d];
565:               }
566:             }
567:           }
568:           cell++;
569:         }
570:       }
571:       for (v = 0; v < numVertices; ++v) {
572:         PetscSectionGetOffset(coordSection, v + trueNumCells, &off);
573:         for (d = 0; d < embedDim; ++d) {
574:           coords[off+d] = coordsIn[periodicMapI[v]*3+d];
575:         }
576:       }
577:     } else {
578:       for (v = 0; v < numVertices; ++v) {
579:         for (d = 0; d < embedDim; ++d) {
580:           coords[v*embedDim+d] = coordsIn[v*3+d];
581:         }
582:       }
583:     }
584:   }
585:   DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);
586:   VecRestoreArray(coordinates, &coords);
587:   PetscFree(coordsIn);
588:   DMSetCoordinatesLocal(*dm, coordinates);
589:   VecDestroy(&coordinates);
590:   PetscFree(periodicMap);
591:   PetscFree(periodicMapI);
592:   PetscBTDestroy(&periodicV);
593:   PetscBTDestroy(&periodicC);
594:   /* Clean up intermediate storage */
595:   if (!rank || binary) PetscFree(gmsh_elem);
596:   PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
597:   return(0);
598: }