Actual source code: plexgmsh.c

petsc-3.6.4 2016-04-12
Report Typos and Errors
  1: #define PETSCDM_DLL
  2: #include <petsc/private/dmpleximpl.h>    /*I   "petscdmplex.h"   I*/

  6: /*@C
  7:   DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file

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

 13:   Output Parameter:
 14: . dm  - The DM object representing the mesh

 16:   Level: beginner

 18: .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
 19: @*/
 20: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
 21: {
 22:   PetscViewer     viewer, vheader;
 23:   PetscMPIInt     rank;
 24:   PetscViewerType vtype;
 25:   char            line[PETSC_MAX_PATH_LEN];
 26:   int             snum;
 27:   PetscBool       match;
 28:   int             fT;
 29:   PetscInt        fileType;
 30:   float           version;
 31:   PetscErrorCode  ierr;

 34:   MPI_Comm_rank(comm, &rank);
 35:   /* Determine Gmsh file type (ASCII or binary) from file header */
 36:   PetscViewerCreate(comm, &vheader);
 37:   PetscViewerSetType(vheader, PETSCVIEWERASCII);
 38:   PetscViewerFileSetMode(vheader, FILE_MODE_READ);
 39:   PetscViewerFileSetName(vheader, filename);
 40:   if (!rank) {
 41:     /* Read only the first two lines of the Gmsh file */
 42:     PetscViewerRead(vheader, line, 1, NULL, PETSC_STRING);
 43:     PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &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, &fT);
 47:     fileType = (PetscInt) fT;
 48:     if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
 49:     if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
 50:   }
 51:   MPI_Bcast(&fileType, 1, MPIU_INT, 0, comm);
 52:   /* Create appropriate viewer and build plex */
 53:   if (fileType == 0) vtype = PETSCVIEWERASCII;
 54:   else vtype = PETSCVIEWERBINARY;
 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:   PetscViewerDestroy(&vheader);
 62:   return(0);
 63: }


 68: /*@
 69:   DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer

 71:   Collective on comm

 73:   Input Parameters:
 74: + comm  - The MPI communicator
 75: . viewer - The Viewer associated with a Gmsh file
 76: - interpolate - Create faces and edges in the mesh

 78:   Output Parameter:
 79: . dm  - The DM object representing the mesh

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

 84:   Level: beginner

 86: .keywords: mesh,Gmsh
 87: .seealso: DMPLEX, DMCreate()
 88: @*/
 89: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
 90: {
 91:   PetscViewerType vtype;
 92:   GmshElement   *gmsh_elem;
 93:   PetscSection   coordSection;
 94:   Vec            coordinates;
 95:   PetscScalar   *coords, *coordsIn = NULL;
 96:   PetscInt       dim = 0, coordSize, c, v, d, cell;
 97:   int            i, numVertices = 0, numCells = 0, trueNumCells = 0, snum;
 98:   PetscMPIInt    num_proc, rank;
 99:   char           line[PETSC_MAX_PATH_LEN];
100:   PetscBool      match, binary, bswap = PETSC_FALSE;

104:   MPI_Comm_rank(comm, &rank);
105:   MPI_Comm_size(comm, &num_proc);
106:   DMCreate(comm, dm);
107:   DMSetType(*dm, DMPLEX);
108:   PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);
109:   PetscViewerGetType(viewer, &vtype);
110:   PetscStrcmp(vtype, PETSCVIEWERBINARY, &binary);
111:   if (!rank || binary) {
112:     PetscBool match;
113:     int       fileType, dataSize;
114:     float     version;

116:     /* Read header */
117:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
118:     PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);
119:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
120:     PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
121:     snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
122:     if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line);
123:     if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0");
124:     if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize);
125:     if (binary) {
126:       int checkInt;
127:       PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);
128:       if (checkInt != 1) {
129:         PetscByteSwap(&checkInt, PETSC_ENUM, 1);
130:         if (checkInt == 1) bswap = PETSC_TRUE;
131:         else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType);
132:       }
133:     } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType);
134:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
135:     PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);
136:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
137:     /* Read vertices */
138:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
139:     PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);
140:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
141:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
142:     snum = sscanf(line, "%d", &numVertices);
143:     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
144:     PetscMalloc1(numVertices*3, &coordsIn);
145:     if (binary) {
146:       size_t doubleSize, intSize;
147:       PetscInt elementSize;
148:       char *buffer;
149:       PetscScalar *baseptr;
150:       PetscDataTypeGetSize(PETSC_ENUM, &intSize);
151:       PetscDataTypeGetSize(PETSC_DOUBLE, &doubleSize);
152:       elementSize = (intSize + 3*doubleSize);
153:       PetscMalloc1(elementSize*numVertices, &buffer);
154:       PetscViewerRead(viewer, buffer, elementSize*numVertices, NULL, PETSC_CHAR);
155:       if (bswap) PetscByteSwap(buffer, PETSC_CHAR, elementSize*numVertices);
156:       for (v = 0; v < numVertices; ++v) {
157:         baseptr = ((PetscScalar*)(buffer+v*elementSize+intSize));
158:         coordsIn[v*3+0] = baseptr[0];
159:         coordsIn[v*3+1] = baseptr[1];
160:         coordsIn[v*3+2] = baseptr[2];
161:       }
162:       PetscFree(buffer);
163:     } else {
164:       for (v = 0; v < numVertices; ++v) {
165:         PetscViewerRead(viewer, &i, 1, NULL, PETSC_ENUM);
166:         PetscViewerRead(viewer, &(coordsIn[v*3]), 3, NULL, PETSC_DOUBLE);
167:         if (i != (int)v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+1);
168:       }
169:     }
170:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
171:     PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);
172:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
173:     /* Read cells */
174:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
175:     PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);
176:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
177:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
178:     snum = sscanf(line, "%d", &numCells);
179:     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
180:   }

182:   if (!rank || binary) {
183:     /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
184:        file contents multiple times to figure out the true number of cells and facets
185:        in the given mesh. To make this more efficient we read the file contents only
186:        once and store them in memory, while determining the true number of cells. */
187:     DMPlexCreateGmsh_ReadElement(viewer, numCells, binary, bswap, &gmsh_elem);
188:     for (trueNumCells=0, c = 0; c < numCells; ++c) {
189:       if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;}
190:       if (gmsh_elem[c].dim == dim) trueNumCells++;
191:     }
192:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);;
193:     PetscStrncmp(line, "$EndElements", PETSC_MAX_PATH_LEN, &match);
194:     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
195:   }
196:   /* For binary we read on all ranks, but only build the plex on rank 0 */
197:   if (binary && rank) {trueNumCells = 0; numVertices = 0;};
198:   /* Allocate the cell-vertex mesh */
199:   DMPlexSetChart(*dm, 0, trueNumCells+numVertices);
200:   if (!rank) {
201:     for (cell = 0, c = 0; c < numCells; ++c) {
202:       if (gmsh_elem[c].dim == dim) {
203:         DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);
204:         cell++;
205:       }
206:     }
207:   }
208:   DMSetUp(*dm);
209:   /* Add cell-vertex connections */
210:   if (!rank) {
211:     PetscInt pcone[8], corner;
212:     for (cell = 0, c = 0; c < numCells; ++c) {
213:       if (gmsh_elem[c].dim == dim) {
214:         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
215:           pcone[corner] = gmsh_elem[c].nodes[corner] + trueNumCells-1;
216:         }
217:         DMPlexSetCone(*dm, cell, pcone);
218:         cell++;
219:       }
220:     }
221:   }
222:   MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
223:   DMSetDimension(*dm, dim);
224:   DMPlexSymmetrize(*dm);
225:   DMPlexStratify(*dm);
226:   if (interpolate) {
227:     DM idm = NULL;

229:     DMPlexInterpolate(*dm, &idm);
230:     DMDestroy(dm);
231:     *dm  = idm;
232:   }

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

238:     DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
239:     for (c = 0; c < numCells; ++c) {
240:       if (gmsh_elem[c].dim == dim-1) {
241:         PetscInt joinSize;
242:         const PetscInt *join;
243:         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
244:           pcone[corner] = gmsh_elem[c].nodes[corner] + vStart - 1;
245:         }
246:         DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
247:         if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id);
248:         DMPlexSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);
249:         DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);
250:       }
251:     }
252:   }

254:   /* Read coordinates */
255:   DMGetCoordinateSection(*dm, &coordSection);
256:   PetscSectionSetNumFields(coordSection, 1);
257:   PetscSectionSetFieldComponents(coordSection, 0, dim);
258:   PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);
259:   for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
260:     PetscSectionSetDof(coordSection, v, dim);
261:     PetscSectionSetFieldDof(coordSection, v, 0, dim);
262:   }
263:   PetscSectionSetUp(coordSection);
264:   PetscSectionGetStorageSize(coordSection, &coordSize);
265:   VecCreate(comm, &coordinates);
266:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
267:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
268:   VecSetType(coordinates, VECSTANDARD);
269:   VecGetArray(coordinates, &coords);
270:   if (!rank) {
271:     for (v = 0; v < numVertices; ++v) {
272:       for (d = 0; d < dim; ++d) {
273:         coords[v*dim+d] = coordsIn[v*3+d];
274:       }
275:     }
276:   }
277:   VecRestoreArray(coordinates, &coords);
278:   PetscFree(coordsIn);
279:   DMSetCoordinatesLocal(*dm, coordinates);
280:   VecDestroy(&coordinates);
281:   /* Clean up intermediate storage */
282:   if (!rank || binary) PetscFree(gmsh_elem);
283:   PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);
284:   return(0);
285: }

289: PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscInt numCells, PetscBool binary, PetscBool byteSwap, GmshElement **gmsh_elems)
290: {
291:   PetscInt       c, p;
292:   GmshElement   *elements;
293:   int            i, cellType, dim, numNodes, numElem, numTags;
294:   int            ibuf[16];

298:   PetscMalloc1(numCells, &elements);
299:   for (c = 0; c < numCells;) {
300:     PetscViewerRead(viewer, &ibuf, 3, NULL, PETSC_ENUM);
301:     if (byteSwap) PetscByteSwap(&ibuf, PETSC_ENUM, 3);
302:     if (binary) {
303:       cellType = ibuf[0];
304:       numElem = ibuf[1];
305:       numTags = ibuf[2];
306:     } else {
307:       elements[c].id = ibuf[0];
308:       cellType = ibuf[1];
309:       numTags = ibuf[2];
310:       numElem = 1;
311:     }
312:     switch (cellType) {
313:     case 1: /* 2-node line */
314:       dim = 1;
315:       numNodes = 2;
316:       break;
317:     case 2: /* 3-node triangle */
318:       dim = 2;
319:       numNodes = 3;
320:       break;
321:     case 3: /* 4-node quadrangle */
322:       dim = 2;
323:       numNodes = 4;
324:       break;
325:     case 4: /* 4-node tetrahedron */
326:       dim  = 3;
327:       numNodes = 4;
328:       break;
329:     case 5: /* 8-node hexahedron */
330:       dim = 3;
331:       numNodes = 8;
332:       break;
333:     case 15: /* 1-node vertex */
334:       dim = 0;
335:       numNodes = 1;
336:       break;
337:     default:
338:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
339:     }
340:     if (binary) {
341:       const PetscInt nint = numNodes + numTags + 1;
342:       for (i = 0; i < numElem; ++i, ++c) {
343:         /* Loop over inner binary element block */
344:         elements[c].dim = dim;
345:         elements[c].numNodes = numNodes;
346:         elements[c].numTags = numTags;

348:         PetscViewerRead(viewer, &ibuf, nint, NULL, PETSC_ENUM);
349:         if (byteSwap) PetscByteSwap( &ibuf, PETSC_ENUM, nint);
350:         elements[c].id = ibuf[0];
351:         for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p];
352:         for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
353:       }
354:     } else {
355:       elements[c].dim = dim;
356:       elements[c].numNodes = numNodes;
357:       elements[c].numTags = numTags;
358:       PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);
359:       PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);
360:       c++;
361:     }
362:   }
363:   *gmsh_elems = elements;
364:   return(0);
365: }