Actual source code: plexgmsh.c

petsc-3.10.5 2019-03-28
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 6: /* 6-node wedge */
128:       dim = 3;
129:       numNodes = 6;
130:       break;
131:     case 8: /* 3-node 2nd order line */
132:       dim = 1;
133:       numNodes = 2;
134:       numNodesIgnore = 1;
135:       break;
136:     case 9: /* 6-node 2nd order triangle */
137:       dim = 2;
138:       numNodes = 3;
139:       numNodesIgnore = 3;
140:       break;
141:     case 13: /* 18-node 2nd wedge */
142:       dim = 3;
143:       numNodes = 6;
144:       numNodesIgnore = 12;
145:       break;
146:     case 15: /* 1-node vertex */
147:       dim = 0;
148:       numNodes = 1;
149:       break;
150:     case 7: /* 5-node pyramid */
151:     case 10: /* 9-node 2nd order quadrangle */
152:     case 11: /* 10-node 2nd order tetrahedron */
153:     case 12: /* 27-node 2nd order hexhedron */
154:     case 14: /* 14-node 2nd order pyramid */
155:     default:
156:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
157:     }
158:     if (binary) {
159:       const int nint = 1 + numTags + numNodes + numNodesIgnore;
160:       /* Loop over element blocks */
161:       for (i = 0; i < numElem; ++i, ++c) {
162:         PetscViewerRead(viewer, ibuf, nint, NULL, PETSC_ENUM);
163:         if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, nint);}
164:         elements[c].dim = dim;
165:         elements[c].numNodes = numNodes;
166:         elements[c].numTags = numTags;
167:         elements[c].id = ibuf[0];
168:         elements[c].cellType = cellType;
169:         for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p];
170:         for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
171:       }
172:     } else {
173:       elements[c].dim = dim;
174:       elements[c].numNodes = numNodes;
175:       elements[c].numTags = numTags;
176:       elements[c].cellType = cellType;
177:       PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);
178:       PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);
179:       PetscViewerRead(viewer, ibuf, numNodesIgnore, NULL, PETSC_ENUM);
180:       c++;
181:     }
182:   }
183:   *gmsh_elems = elements;
184:   return(0);
185: }

187: /*@
188:   DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer

190:   Collective on comm

192:   Input Parameters:
193: + comm  - The MPI communicator
194: . viewer - The Viewer associated with a Gmsh file
195: - interpolate - Create faces and edges in the mesh

197:   Output Parameter:
198: . dm  - The DM object representing the mesh

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

203:   Level: beginner

205: .keywords: mesh,Gmsh
206: .seealso: DMPLEX, DMCreate()
207: @*/
208: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
209: {
210:   PetscViewer    parentviewer = NULL;
211:   double        *coordsIn = NULL;
212:   GmshElement   *gmsh_elem = NULL;
213:   PetscSection   coordSection;
214:   Vec            coordinates;
215:   PetscBT        periodicV = NULL, periodicC = NULL;
216:   PetscScalar   *coords;
217:   PetscInt       dim = 0, embedDim = -1, coordSize, c, v, d, r, cell, *periodicMap = NULL, *periodicMapI = NULL, *hybridMap = NULL, cMax = PETSC_DETERMINE;
218:   int            i, numVertices = 0, numCells = 0, trueNumCells = 0, numRegions = 0, snum, shift = 1;
219:   PetscMPIInt    rank;
220:   char           line[PETSC_MAX_PATH_LEN];
221:   PetscBool      binary, byteSwap = PETSC_FALSE, zerobase = PETSC_FALSE, periodic = PETSC_FALSE, usemarker = PETSC_FALSE;
222:   PetscBool      enable_hybrid = PETSC_FALSE;

226:   MPI_Comm_rank(comm, &rank);
227:   PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_hybrid", &enable_hybrid, NULL);
228:   PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_periodic", &periodic, NULL);
229:   PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_use_marker", &usemarker, NULL);
230:   PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_zero_base", &zerobase, NULL);
231:   PetscOptionsGetInt(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_spacedim", &embedDim, NULL);
232:   if (zerobase) shift = 0;

234:   DMCreate(comm, dm);
235:   DMSetType(*dm, DMPLEX);
236:   PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);

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

240:   /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
241:   if (binary) {
242:     parentviewer = viewer;
243:     PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
244:   }

246:   if (!rank) {
247:     PetscBool match, hybrid;
248:     int       fileType, dataSize;
249:     float     version;

251:     /* Read header */
252:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
253:     PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);
254:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
255:     PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
256:     snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
257:     if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
258:     if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
259:     if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
260:     if (binary) {
261:       int checkInt;
262:       PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);
263:       if (checkInt != 1) {
264:         PetscByteSwap(&checkInt, PETSC_ENUM, 1);
265:         if (checkInt == 1) byteSwap = PETSC_TRUE;
266:         else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType);
267:       }
268:     } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType);
269:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
270:     PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);
271:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");

273:     /* OPTIONAL Read physical names */
274:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
275:     PetscStrncmp(line, "$PhysicalNames", PETSC_MAX_PATH_LEN, &match);
276:     if (match) {
277:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
278:       snum = sscanf(line, "%d", &numRegions);
279:       if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
280:       for (r = 0; r < numRegions; ++r) {
281:         PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
282:       }
283:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
284:       PetscStrncmp(line, "$EndPhysicalNames", PETSC_MAX_PATH_LEN, &match);
285:       if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
286:       /* Initial read for vertex section */
287:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
288:     }

290:     /* Read vertices */
291:     PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);
292:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
293:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
294:     snum = sscanf(line, "%d", &numVertices);
295:     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
296:     DMPlexCreateGmsh_ReadNodes(viewer, binary, byteSwap, shift, numVertices, &coordsIn);
297:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
298:     PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);
299:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");

301:     /* Read cells */
302:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
303:     PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);
304:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
305:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
306:     snum = sscanf(line, "%d", &numCells);
307:     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
308:     /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
309:        file contents multiple times to figure out the true number of cells and facets
310:        in the given mesh. To make this more efficient we read the file contents only
311:        once and store them in memory, while determining the true number of cells. */
312:     DMPlexCreateGmsh_ReadElements(viewer, binary, byteSwap, shift, numCells, &gmsh_elem);
313:     hybrid = PETSC_FALSE;
314:     for (trueNumCells = 0, c = 0; c < numCells; ++c) {
315:       int on = -1;
316:       if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;}
317:       if (gmsh_elem[c].dim == dim) {hybrid = (trueNumCells ? (on != gmsh_elem[c].cellType ? on = gmsh_elem[c].cellType,PETSC_TRUE : hybrid) : (on = gmsh_elem[c].cellType, PETSC_FALSE) ); trueNumCells++;}
318:       /* wedges always indicate an hybrid mesh in PLEX */
319:       if (on == 6 || on == 13) hybrid = PETSC_TRUE;
320:     }
321:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
322:     PetscStrncmp(line, "$EndElements", 12, &match);
323:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");

325:     /* Renumber cells for hybrid grids */
326:     if (hybrid && enable_hybrid) {
327:       PetscInt hc1 = 0, hc2 = 0, *hybridCells1 = NULL, *hybridCells2 = NULL;
328:       PetscInt cell, tn, *tp;
329:       int      n1 = 0,n2 = 0;

331:       PetscMalloc1(trueNumCells, &hybridCells1);
332:       PetscMalloc1(trueNumCells, &hybridCells2);
333:       for (cell = 0, c = 0; c < numCells; ++c) {
334:         if (gmsh_elem[c].dim == dim) {
335:           if (!n1) n1 = gmsh_elem[c].cellType;
336:           else if (!n2 && n1 != gmsh_elem[c].cellType) n2 = gmsh_elem[c].cellType;

338:           if      (gmsh_elem[c].cellType == n1) hybridCells1[hc1++] = cell;
339:           else if (gmsh_elem[c].cellType == n2) hybridCells2[hc2++] = cell;
340:           else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle more than 2 cell types");
341:           cell++;
342:         }
343:       }

345:       switch (n1) {
346:       case 2: /* triangles */
347:       case 9:
348:         switch (n2) {
349:         case 0: /* single type mesh */
350:         case 3: /* quads */
351:         case 10:
352:           break;
353:         default:
354:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
355:         }
356:         break;
357:       case 3:
358:       case 10:
359:         switch (n2) {
360:         case 0: /* single type mesh */
361:         case 2: /* swap since we list simplices first */
362:         case 9:
363:           tn  = hc1;
364:           hc1 = hc2;
365:           hc2 = tn;

367:           tp           = hybridCells1;
368:           hybridCells1 = hybridCells2;
369:           hybridCells2 = tp;
370:           break;
371:         default:
372:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
373:         }
374:         break;
375:       case 4: /* tetrahedra */
376:       case 11:
377:         switch (n2) {
378:         case 0: /* single type mesh */
379:         case 6: /* wedges */
380:         case 13:
381:           break;
382:         default:
383:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
384:         }
385:         break;
386:       case 6:
387:       case 13:
388:         switch (n2) {
389:         case 0: /* single type mesh */
390:         case 4: /* swap since we list simplices first */
391:         case 11:
392:           tn  = hc1;
393:           hc1 = hc2;
394:           hc2 = tn;

396:           tp           = hybridCells1;
397:           hybridCells1 = hybridCells2;
398:           hybridCells2 = tp;
399:           break;
400:         default:
401:           SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
402:         }
403:         break;
404:       default:
405:         SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2);
406:       }
407:       cMax = hc1;
408:       PetscMalloc1(trueNumCells, &hybridMap);
409:       for (cell = 0; cell < hc1; cell++) hybridMap[hybridCells1[cell]] = cell;
410:       for (cell = 0; cell < hc2; cell++) hybridMap[hybridCells2[cell]] = cell + hc1;
411:       PetscFree(hybridCells1);
412:       PetscFree(hybridCells2);
413:     }

415:     /* OPTIONAL Read periodic section */
416:     if (periodic) {
417:       PetscInt pVert, *periodicMapT, *aux;
418:       int      numPeriodic;

420:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
421:       PetscStrncmp(line, "$Periodic", 9, &match);
422:       if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
423:       PetscMalloc1(numVertices, &periodicMapT);
424:       PetscBTCreate(numVertices, &periodicV);
425:       for (i = 0; i < numVertices; i++) periodicMapT[i] = i;
426:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
427:       snum = sscanf(line, "%d", &numPeriodic);
428:       if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
429:       for (i = 0; i < numPeriodic; i++) {
430:         int j, edim = -1, slaveTag = -1, masterTag = -1, nNodes, slaveNode, masterNode;

432:         PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
433:         snum = sscanf(line, "%d %d %d", &edim, &slaveTag, &masterTag); /* slaveTag and masterTag are unused */
434:         if (snum != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
435:         PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
436:         snum = sscanf(line, "%d", &nNodes);
437:         if (snum != 1) { /* discard tranformation and try again */
438:           PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
439:           PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
440:           snum = sscanf(line, "%d", &nNodes);
441:           if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
442:         }
443:         for (j = 0; j < nNodes; j++) {
444:           PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
445:           snum = sscanf(line, "%d %d", &slaveNode, &masterNode);
446:           if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
447:           periodicMapT[slaveNode - shift] = masterNode - shift;
448:           PetscBTSet(periodicV, slaveNode - shift);
449:           PetscBTSet(periodicV, masterNode - shift);
450:         }
451:       }
452:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
453:       PetscStrncmp(line, "$EndPeriodic", 12, &match);
454:       if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
455:       /* we may have slaves of slaves */
456:       for (i = 0; i < numVertices; i++) {
457:         while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) {
458:           periodicMapT[i] = periodicMapT[periodicMapT[i]];
459:         }
460:       }
461:       /* periodicMap : from old to new numbering (periodic vertices excluded)
462:          periodicMapI: from new to old numbering */
463:       PetscMalloc1(numVertices, &periodicMap);
464:       PetscMalloc1(numVertices, &periodicMapI);
465:       PetscMalloc1(numVertices, &aux);
466:       for (i = 0, pVert = 0; i < numVertices; i++) {
467:         if (periodicMapT[i] != i) {
468:           pVert++;
469:         } else {
470:           aux[i] = i - pVert;
471:           periodicMapI[i - pVert] = i;
472:         }
473:       }
474:       for (i = 0 ; i < numVertices; i++) {
475:         periodicMap[i] = aux[periodicMapT[i]];
476:       }
477:       PetscFree(periodicMapT);
478:       PetscFree(aux);
479:       /* remove periodic vertices */
480:       numVertices = numVertices - pVert;
481:     }
482:   }

484:   if (parentviewer) {
485:     PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
486:   }

488:   /* Allocate the cell-vertex mesh */
489:   DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
490:   for (cell = 0, c = 0; c < numCells; ++c) {
491:     if (gmsh_elem[c].dim == dim) {
492:       DMPlexSetConeSize(*dm, hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].numNodes);
493:       cell++;
494:     }
495:   }
496:   DMSetUp(*dm);
497:   /* Add cell-vertex connections */
498:   for (cell = 0, c = 0; c < numCells; ++c) {
499:     if (gmsh_elem[c].dim == dim) {
500:       PetscInt pcone[8], corner;
501:       for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
502:         const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
503:         pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells;
504:       }
505:       if (dim == 3) {
506:         /* Tetrahedra are inverted */
507:         if (gmsh_elem[c].cellType == 4) {
508:           PetscInt tmp = pcone[0];
509:           pcone[0] = pcone[1];
510:           pcone[1] = tmp;
511:         }
512:         /* Hexahedra are inverted */
513:         if (gmsh_elem[c].cellType == 5) {
514:           PetscInt tmp = pcone[1];
515:           pcone[1] = pcone[3];
516:           pcone[3] = tmp;
517:         }
518:         /* Prisms are inverted */
519:         if (gmsh_elem[c].cellType == 6) {
520:           PetscInt tmp;

522:           tmp      = pcone[1];
523:           pcone[1] = pcone[2];
524:           pcone[2] = tmp;
525:           tmp      = pcone[4];
526:           pcone[4] = pcone[5];
527:           pcone[5] = tmp;
528:         }
529:       } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */
530:         /* quads are input to PLEX as prisms */
531:         if (gmsh_elem[c].cellType == 3) {
532:           PetscInt tmp = pcone[2];
533:           pcone[2] = pcone[3];
534:           pcone[3] = tmp;
535:         }
536:       }
537:       DMPlexSetCone(*dm, hybridMap ? hybridMap[cell] : cell, pcone);
538:       cell++;
539:     }
540:   }
541:   MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
542:   DMSetDimension(*dm, dim);
543:   DMPlexSetHybridBounds(*dm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);
544:   DMPlexSymmetrize(*dm);
545:   DMPlexStratify(*dm);
546:   if (interpolate) {
547:     DM idm;

549:     DMPlexInterpolate(*dm, &idm);
550:     DMDestroy(dm);
551:     *dm  = idm;
552:   }

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

558:     DMCreateLabel(*dm, "marker");
559:     DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
560:     for (f = fStart; f < fEnd; ++f) {
561:       PetscInt suppSize;

563:       DMPlexGetSupportSize(*dm, f, &suppSize);
564:       if (suppSize == 1) {
565:         PetscInt *cone = NULL, coneSize, p;

567:         DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
568:         for (p = 0; p < coneSize; p += 2) {
569:           DMSetLabelValue(*dm, "marker", cone[p], 1);
570:         }
571:         DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
572:       }
573:     }
574:   }

576:   if (!rank) {
577:     PetscInt vStart, vEnd;

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

582:       /* Create face sets */
583:       if (interpolate && gmsh_elem[c].dim == dim-1) {
584:         const PetscInt *join;
585:         PetscInt        joinSize, pcone[8], corner;
586:         /* Find the relevant facet with vertex joins */
587:         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
588:           const PetscInt cc = gmsh_elem[c].nodes[corner] - shift;
589:           pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart;
590:         }
591:         DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, pcone, &joinSize, &join);
592:         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);
593:         DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
594:         DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
595:       }

597:       /* Create cell sets */
598:       if (gmsh_elem[c].dim == dim) {
599:         if (gmsh_elem[c].numTags > 0) {
600:           DMSetLabelValue(*dm, "Cell Sets", hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].tags[0]);
601:         }
602:         cell++;
603:       }

605:       /* Create vertex sets */
606:       if (gmsh_elem[c].dim == 0) {
607:         if (gmsh_elem[c].numTags > 0) {
608:           const PetscInt cc = gmsh_elem[c].nodes[0] - shift;
609:           const PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart;
610:           DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);
611:         }
612:       }
613:     }
614:   }

616:   /* Create coordinates */
617:   if (embedDim < 0) embedDim = dim;
618:   DMSetCoordinateDim(*dm, embedDim);
619:   DMGetCoordinateSection(*dm, &coordSection);
620:   PetscSectionSetNumFields(coordSection, 1);
621:   PetscSectionSetFieldComponents(coordSection, 0, embedDim);
622:   if (periodicMap) { /* we need to localize coordinates on cells */
623:     PetscSectionSetChart(coordSection, 0, trueNumCells + numVertices);
624:   } else {
625:     PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);
626:   }
627:   for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
628:     PetscSectionSetDof(coordSection, v, embedDim);
629:     PetscSectionSetFieldDof(coordSection, v, 0, embedDim);
630:   }
631:   if (periodicMap) {
632:     PetscBTCreate(trueNumCells, &periodicC);
633:     for (cell = 0, c = 0; c < numCells; ++c) {
634:       if (gmsh_elem[c].dim == dim) {
635:         PetscInt  corner;
636:         PetscBool pc = PETSC_FALSE;
637:         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
638:           pc = (PetscBool)(pc || PetscBTLookup(periodicV, gmsh_elem[c].nodes[corner] - shift));
639:         }
640:         if (pc) {
641:           PetscInt dof = gmsh_elem[c].numNodes*embedDim;
642:           PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
643:           PetscBTSet(periodicC, ucell);
644:           PetscSectionSetDof(coordSection, ucell, dof);
645:           PetscSectionSetFieldDof(coordSection, ucell, 0, dof);
646:         }
647:         cell++;
648:       }
649:     }
650:   }
651:   PetscSectionSetUp(coordSection);
652:   PetscSectionGetStorageSize(coordSection, &coordSize);
653:   VecCreate(PETSC_COMM_SELF, &coordinates);
654:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
655:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
656:   VecSetBlockSize(coordinates, embedDim);
657:   VecSetType(coordinates, VECSTANDARD);
658:   VecGetArray(coordinates, &coords);
659:   if (periodicMap) {
660:     PetscInt off;

662:     for (cell = 0, c = 0; c < numCells; ++c) {
663:       PetscInt pcone[8], corner;
664:       if (gmsh_elem[c].dim == dim) {
665:         PetscInt ucell = hybridMap ? hybridMap[cell] : cell;
666:         if (PetscUnlikely(PetscBTLookup(periodicC, ucell))) {
667:           for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
668:             pcone[corner] = gmsh_elem[c].nodes[corner] - shift;
669:           }
670:           if (dim == 3) {
671:             /* Tetrahedra are inverted */
672:             if (gmsh_elem[c].cellType == 4) {
673:               PetscInt tmp = pcone[0];
674:               pcone[0] = pcone[1];
675:               pcone[1] = tmp;
676:             }
677:             /* Hexahedra are inverted */
678:             if (gmsh_elem[c].cellType == 5) {
679:               PetscInt tmp = pcone[1];
680:               pcone[1] = pcone[3];
681:               pcone[3] = tmp;
682:             }
683:             /* Prisms are inverted */
684:             if (gmsh_elem[c].cellType == 6) {
685:               PetscInt tmp;

687:               tmp      = pcone[1];
688:               pcone[1] = pcone[2];
689:               pcone[2] = tmp;
690:               tmp      = pcone[4];
691:               pcone[4] = pcone[5];
692:               pcone[5] = tmp;
693:             }
694:           } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */
695:             /* quads are input to PLEX as prisms */
696:             if (gmsh_elem[c].cellType == 3) {
697:               PetscInt tmp = pcone[2];
698:               pcone[2] = pcone[3];
699:               pcone[3] = tmp;
700:             }
701:           }
702:           PetscSectionGetOffset(coordSection, ucell, &off);
703:           for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
704:             v = pcone[corner];
705:             for (d = 0; d < embedDim; ++d) {
706:               coords[off++] = (PetscReal) coordsIn[v*3+d];
707:             }
708:           }
709:         }
710:         cell++;
711:       }
712:     }
713:     for (v = 0; v < numVertices; ++v) {
714:       PetscSectionGetOffset(coordSection, v + trueNumCells, &off);
715:       for (d = 0; d < embedDim; ++d) {
716:         coords[off+d] = (PetscReal) coordsIn[periodicMapI[v]*3+d];
717:       }
718:     }
719:   } else {
720:     for (v = 0; v < numVertices; ++v) {
721:       for (d = 0; d < embedDim; ++d) {
722:         coords[v*embedDim+d] = (PetscReal) coordsIn[v*3+d];
723:       }
724:     }
725:   }
726:   VecRestoreArray(coordinates, &coords);
727:   DMSetCoordinatesLocal(*dm, coordinates);
728:   DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);

730:   PetscFree(coordsIn);
731:   PetscFree(gmsh_elem);
732:   PetscFree(hybridMap);
733:   PetscFree(periodicMap);
734:   PetscFree(periodicMapI);
735:   PetscBTDestroy(&periodicV);
736:   PetscBTDestroy(&periodicC);
737:   VecDestroy(&coordinates);

739:   PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
740:   return(0);
741: }