xref: /petsc/src/dm/impls/plex/plexgmsh.c (revision 0c070f12aa33cd544b9d7968a9e5c044728f8926)
1 #define PETSCDM_DLL
2 #include <petsc/private/dmpleximpl.h>    /*I   "petscdmplex.h"   I*/
3 
4 #undef __FUNCT__
5 #define __FUNCT__ "DMPlexCreateGmshFromFile"
6 /*@C
7   DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file
8 
9 + comm        - The MPI communicator
10 . filename    - Name of the Gmsh file
11 - interpolate - Create faces and edges in the mesh
12 
13   Output Parameter:
14 . dm  - The DM object representing the mesh
15 
16   Level: beginner
17 
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;
32 
33   PetscFunctionBegin;
34   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
35   /* Determine Gmsh file type (ASCII or binary) from file header */
36   ierr = PetscViewerCreate(comm, &vheader);CHKERRQ(ierr);
37   ierr = PetscViewerSetType(vheader, PETSCVIEWERASCII);CHKERRQ(ierr);
38   ierr = PetscViewerFileSetMode(vheader, FILE_MODE_READ);CHKERRQ(ierr);
39   ierr = PetscViewerFileSetName(vheader, filename);CHKERRQ(ierr);
40   if (!rank) {
41     /* Read only the first two lines of the Gmsh file */
42     ierr = PetscViewerRead(vheader, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
43     ierr = PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
44     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
45     ierr = PetscViewerRead(vheader, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr);
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   ierr = MPI_Bcast(&fileType, 1, MPIU_INT, 0, comm);CHKERRQ(ierr);
52   /* Create appropriate viewer and build plex */
53   if (fileType == 0) vtype = PETSCVIEWERASCII;
54   else vtype = PETSCVIEWERBINARY;
55   ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr);
56   ierr = PetscViewerSetType(viewer, vtype);CHKERRQ(ierr);
57   ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr);
58   ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr);
59   ierr = DMPlexCreateGmsh(comm, viewer, interpolate, dm);CHKERRQ(ierr);
60   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
61   ierr = PetscViewerDestroy(&vheader);CHKERRQ(ierr);
62   PetscFunctionReturn(0);
63 }
64 
65 
66 #undef __FUNCT__
67 #define __FUNCT__ "DMPlexCreateGmsh"
68 /*@
69   DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer
70 
71   Collective on comm
72 
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
77 
78   Output Parameter:
79 . dm  - The DM object representing the mesh
80 
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
83 
84   Level: beginner
85 
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, r, cell;
97   int            i, numVertices = 0, numCells = 0, trueNumCells = 0, numRegions = 0, snum;
98   PetscMPIInt    num_proc, rank;
99   char           line[PETSC_MAX_PATH_LEN];
100   PetscBool      match, binary, bswap = PETSC_FALSE;
101   PetscErrorCode ierr;
102 
103   PetscFunctionBegin;
104   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
105   ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr);
106   ierr = DMCreate(comm, dm);CHKERRQ(ierr);
107   ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
108   ierr = PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr);
109   ierr = PetscViewerGetType(viewer, &vtype);CHKERRQ(ierr);
110   ierr = PetscStrcmp(vtype, PETSCVIEWERBINARY, &binary);CHKERRQ(ierr);
111   if (!rank || binary) {
112     PetscBool match;
113     int       fileType, dataSize;
114     float     version;
115 
116     /* Read header */
117     ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
118     ierr = PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
119     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
120     ierr = PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);CHKERRQ(ierr);
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       ierr = PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);CHKERRQ(ierr);
128       if (checkInt != 1) {
129         ierr = PetscByteSwap(&checkInt, PETSC_ENUM, 1);CHKERRQ(ierr);
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     ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
135     ierr = PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
136     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
137     /* OPTIONAL Read physical names */
138     ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
139     ierr = PetscStrncmp(line, "$PhysicalNames", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
140     if (match) {
141       ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
142       snum = sscanf(line, "%d", &numRegions);
143       if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
144       for (r = 0; r < numRegions; ++r) {
145         PetscInt rdim, tag;
146 
147         ierr = PetscViewerRead(viewer, &rdim, 1, NULL, PETSC_ENUM);CHKERRQ(ierr);
148         ierr = PetscViewerRead(viewer, &tag,  1, NULL, PETSC_ENUM);CHKERRQ(ierr);
149         ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
150       }
151       ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
152       ierr = PetscStrncmp(line, "$EndPhysicalNames", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
153       if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
154       /* Initial read for vertex section */
155       ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
156     }
157     /* Read vertices */
158     ierr = PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
159     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
160     ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
161     snum = sscanf(line, "%d", &numVertices);
162     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
163     ierr = PetscMalloc1(numVertices*3, &coordsIn);CHKERRQ(ierr);
164     if (binary) {
165       size_t doubleSize, intSize;
166       PetscInt elementSize;
167       char *buffer;
168       PetscScalar *baseptr;
169       ierr = PetscDataTypeGetSize(PETSC_ENUM, &intSize);CHKERRQ(ierr);
170       ierr = PetscDataTypeGetSize(PETSC_DOUBLE, &doubleSize);CHKERRQ(ierr);
171       elementSize = (intSize + 3*doubleSize);
172       ierr = PetscMalloc1(elementSize*numVertices, &buffer);CHKERRQ(ierr);
173       ierr = PetscViewerRead(viewer, buffer, elementSize*numVertices, NULL, PETSC_CHAR);CHKERRQ(ierr);
174       if (bswap) ierr = PetscByteSwap(buffer, PETSC_CHAR, elementSize*numVertices);CHKERRQ(ierr);
175       for (v = 0; v < numVertices; ++v) {
176         baseptr = ((PetscScalar*)(buffer+v*elementSize+intSize));
177         coordsIn[v*3+0] = baseptr[0];
178         coordsIn[v*3+1] = baseptr[1];
179         coordsIn[v*3+2] = baseptr[2];
180       }
181       ierr = PetscFree(buffer);CHKERRQ(ierr);
182     } else {
183       for (v = 0; v < numVertices; ++v) {
184         ierr = PetscViewerRead(viewer, &i, 1, NULL, PETSC_ENUM);CHKERRQ(ierr);
185         ierr = PetscViewerRead(viewer, &(coordsIn[v*3]), 3, NULL, PETSC_DOUBLE);CHKERRQ(ierr);
186         if (i != (int)v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+1);
187       }
188     }
189     ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);;
190     ierr = PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
191     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
192     /* Read cells */
193     ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);;
194     ierr = PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
195     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
196     ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);;
197     snum = sscanf(line, "%d", &numCells);
198     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
199   }
200 
201   if (!rank || binary) {
202     /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
203        file contents multiple times to figure out the true number of cells and facets
204        in the given mesh. To make this more efficient we read the file contents only
205        once and store them in memory, while determining the true number of cells. */
206     ierr = DMPlexCreateGmsh_ReadElement(viewer, numCells, binary, bswap, &gmsh_elem);CHKERRQ(ierr);
207     for (trueNumCells=0, c = 0; c < numCells; ++c) {
208       if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;}
209       if (gmsh_elem[c].dim == dim) trueNumCells++;
210     }
211     ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);;
212     ierr = PetscStrncmp(line, "$EndElements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
213     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
214   }
215   /* For binary we read on all ranks, but only build the plex on rank 0 */
216   if (binary && rank) {trueNumCells = 0; numVertices = 0;};
217   /* Allocate the cell-vertex mesh */
218   ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr);
219   if (!rank) {
220     for (cell = 0, c = 0; c < numCells; ++c) {
221       if (gmsh_elem[c].dim == dim) {
222         ierr = DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);CHKERRQ(ierr);
223         cell++;
224       }
225     }
226   }
227   ierr = DMSetUp(*dm);CHKERRQ(ierr);
228   /* Add cell-vertex connections */
229   if (!rank) {
230     PetscInt pcone[8], corner;
231     for (cell = 0, c = 0; c < numCells; ++c) {
232       if (gmsh_elem[c].dim == dim) {
233         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
234           pcone[corner] = gmsh_elem[c].nodes[corner] + trueNumCells-1;
235         }
236         if (dim == 3) {
237           /* Tetrahedra are inverted */
238           if (gmsh_elem[c].numNodes == 4) {
239             PetscInt tmp = pcone[0];
240             pcone[0] = pcone[1];
241             pcone[1] = tmp;
242           }
243           /* Hexahedra are inverted */
244           if (gmsh_elem[c].numNodes == 8) {
245             PetscInt tmp = pcone[1];
246             pcone[1] = pcone[3];
247             pcone[3] = tmp;
248           }
249         }
250         ierr = DMPlexSetCone(*dm, cell, pcone);CHKERRQ(ierr);
251         cell++;
252       }
253     }
254   }
255   ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr);
256   ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
257   ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
258   ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
259   if (interpolate) {
260     DM idm = NULL;
261 
262     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
263     ierr = DMDestroy(dm);CHKERRQ(ierr);
264     *dm  = idm;
265   }
266 
267   if (!rank) {
268     /* Apply boundary IDs by finding the relevant facets with vertex joins */
269     PetscInt pcone[8], corner, vStart, vEnd;
270 
271     ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
272     for (c = 0; c < numCells; ++c) {
273       if (gmsh_elem[c].dim == dim-1) {
274         PetscInt joinSize;
275         const PetscInt *join;
276         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
277           pcone[corner] = gmsh_elem[c].nodes[corner] + vStart - 1;
278         }
279         ierr = DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr);
280         if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id);
281         ierr = DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);CHKERRQ(ierr);
282         ierr = DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr);
283       }
284     }
285 
286     /* Create cell sets */
287     for (cell = 0, c = 0; c < numCells; ++c) {
288       if (gmsh_elem[c].dim == dim) {
289         if (gmsh_elem[c].numTags > 0) {
290           ierr = DMSetLabelValue(*dm, "Cell Sets", cell, gmsh_elem[c].tags[0]);CHKERRQ(ierr);
291           cell++;
292         }
293       }
294     }
295 
296     /* Create vertex sets */
297     for (c = 0; c < numCells; ++c) {
298       if (gmsh_elem[c].dim == 0) {
299         if (gmsh_elem[c].numTags > 0) {
300           ierr = DMSetLabelValue(*dm, "Vertex Sets", gmsh_elem[c].nodes[0] + vStart - 1, gmsh_elem[c].tags[0]);CHKERRQ(ierr);
301         }
302       }
303     }
304   }
305 
306   /* Read coordinates */
307   ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
308   ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
309   ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr);
310   ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr);
311   for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
312     ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr);
313     ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr);
314   }
315   ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
316   ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
317   ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
318   ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
319   ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
320   ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr);
321   ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr);
322   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
323   if (!rank) {
324     for (v = 0; v < numVertices; ++v) {
325       for (d = 0; d < dim; ++d) {
326         coords[v*dim+d] = coordsIn[v*3+d];
327       }
328     }
329   }
330   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
331   ierr = PetscFree(coordsIn);CHKERRQ(ierr);
332   ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr);
333   ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
334   /* Clean up intermediate storage */
335   if (!rank || binary) ierr = PetscFree(gmsh_elem);CHKERRQ(ierr);
336   ierr = PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr);
337   PetscFunctionReturn(0);
338 }
339 
340 #undef __FUNCT__
341 #define __FUNCT__ "DMPlexCreateGmsh_ReadElement"
342 PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscInt numCells, PetscBool binary, PetscBool byteSwap, GmshElement **gmsh_elems)
343 {
344   PetscInt       c, p;
345   GmshElement   *elements;
346   int            i, cellType, dim, numNodes, numElem, numTags;
347   int            ibuf[16];
348   PetscErrorCode ierr;
349 
350   PetscFunctionBegin;
351   ierr = PetscMalloc1(numCells, &elements);CHKERRQ(ierr);
352   for (c = 0; c < numCells;) {
353     ierr = PetscViewerRead(viewer, &ibuf, 3, NULL, PETSC_ENUM);CHKERRQ(ierr);
354     if (byteSwap) ierr = PetscByteSwap(&ibuf, PETSC_ENUM, 3);CHKERRQ(ierr);
355     if (binary) {
356       cellType = ibuf[0];
357       numElem = ibuf[1];
358       numTags = ibuf[2];
359     } else {
360       elements[c].id = ibuf[0];
361       cellType = ibuf[1];
362       numTags = ibuf[2];
363       numElem = 1;
364     }
365     switch (cellType) {
366     case 1: /* 2-node line */
367       dim = 1;
368       numNodes = 2;
369       break;
370     case 2: /* 3-node triangle */
371       dim = 2;
372       numNodes = 3;
373       break;
374     case 3: /* 4-node quadrangle */
375       dim = 2;
376       numNodes = 4;
377       break;
378     case 4: /* 4-node tetrahedron */
379       dim  = 3;
380       numNodes = 4;
381       break;
382     case 5: /* 8-node hexahedron */
383       dim = 3;
384       numNodes = 8;
385       break;
386     case 15: /* 1-node vertex */
387       dim = 0;
388       numNodes = 1;
389       break;
390     default:
391       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
392     }
393     if (binary) {
394       const PetscInt nint = numNodes + numTags + 1;
395       for (i = 0; i < numElem; ++i, ++c) {
396         /* Loop over inner binary element block */
397         elements[c].dim = dim;
398         elements[c].numNodes = numNodes;
399         elements[c].numTags = numTags;
400 
401         ierr = PetscViewerRead(viewer, &ibuf, nint, NULL, PETSC_ENUM);CHKERRQ(ierr);
402         if (byteSwap) ierr = PetscByteSwap( &ibuf, PETSC_ENUM, nint);CHKERRQ(ierr);
403         elements[c].id = ibuf[0];
404         for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p];
405         for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
406       }
407     } else {
408       elements[c].dim = dim;
409       elements[c].numNodes = numNodes;
410       elements[c].numTags = numTags;
411       ierr = PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);CHKERRQ(ierr);
412       ierr = PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);CHKERRQ(ierr);
413       c++;
414     }
415   }
416   *gmsh_elems = elements;
417   PetscFunctionReturn(0);
418 }
419