xref: /petsc/src/dm/impls/plex/plexgmsh.c (revision ba1e01c44f2f740c97c9026fe63e360558b7709b)
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) {ierr = PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);CHKERRQ(ierr);}
145       ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
146       ierr = PetscStrncmp(line, "$EndPhysicalNames", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
147       if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
148       /* Initial read for vertex section */
149       ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
150     }
151     /* Read vertices */
152     ierr = PetscStrncmp(line, "$Nodes", 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     ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
155     snum = sscanf(line, "%d", &numVertices);
156     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
157     ierr = PetscMalloc1(numVertices*3, &coordsIn);CHKERRQ(ierr);
158     if (binary) {
159       size_t doubleSize, intSize;
160       PetscInt elementSize;
161       char *buffer;
162       PetscScalar *baseptr;
163       ierr = PetscDataTypeGetSize(PETSC_ENUM, &intSize);CHKERRQ(ierr);
164       ierr = PetscDataTypeGetSize(PETSC_DOUBLE, &doubleSize);CHKERRQ(ierr);
165       elementSize = (intSize + 3*doubleSize);
166       ierr = PetscMalloc1(elementSize*numVertices, &buffer);CHKERRQ(ierr);
167       ierr = PetscViewerRead(viewer, buffer, elementSize*numVertices, NULL, PETSC_CHAR);CHKERRQ(ierr);
168       if (bswap) ierr = PetscByteSwap(buffer, PETSC_CHAR, elementSize*numVertices);CHKERRQ(ierr);
169       for (v = 0; v < numVertices; ++v) {
170         baseptr = ((PetscScalar*)(buffer+v*elementSize+intSize));
171         coordsIn[v*3+0] = baseptr[0];
172         coordsIn[v*3+1] = baseptr[1];
173         coordsIn[v*3+2] = baseptr[2];
174       }
175       ierr = PetscFree(buffer);CHKERRQ(ierr);
176     } else {
177       for (v = 0; v < numVertices; ++v) {
178         ierr = PetscViewerRead(viewer, &i, 1, NULL, PETSC_ENUM);CHKERRQ(ierr);
179         ierr = PetscViewerRead(viewer, &(coordsIn[v*3]), 3, NULL, PETSC_DOUBLE);CHKERRQ(ierr);
180         if (i != (int)v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %d should be %d", i, v+1);
181       }
182     }
183     ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);;
184     ierr = PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
185     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
186     /* Read cells */
187     ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);;
188     ierr = PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
189     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
190     ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);;
191     snum = sscanf(line, "%d", &numCells);
192     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
193   }
194 
195   if (!rank || binary) {
196     /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
197        file contents multiple times to figure out the true number of cells and facets
198        in the given mesh. To make this more efficient we read the file contents only
199        once and store them in memory, while determining the true number of cells. */
200     ierr = DMPlexCreateGmsh_ReadElement(viewer, numCells, binary, bswap, &gmsh_elem);CHKERRQ(ierr);
201     for (trueNumCells=0, c = 0; c < numCells; ++c) {
202       if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;}
203       if (gmsh_elem[c].dim == dim) trueNumCells++;
204     }
205     ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);;
206     ierr = PetscStrncmp(line, "$EndElements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr);
207     if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file");
208   }
209   /* For binary we read on all ranks, but only build the plex on rank 0 */
210   if (binary && rank) {trueNumCells = 0; numVertices = 0;};
211   /* Allocate the cell-vertex mesh */
212   ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr);
213   if (!rank) {
214     for (cell = 0, c = 0; c < numCells; ++c) {
215       if (gmsh_elem[c].dim == dim) {
216         ierr = DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);CHKERRQ(ierr);
217         cell++;
218       }
219     }
220   }
221   ierr = DMSetUp(*dm);CHKERRQ(ierr);
222   /* Add cell-vertex connections */
223   if (!rank) {
224     PetscInt pcone[8], corner;
225     for (cell = 0, c = 0; c < numCells; ++c) {
226       if (gmsh_elem[c].dim == dim) {
227         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
228           pcone[corner] = gmsh_elem[c].nodes[corner] + trueNumCells-1;
229         }
230         if (dim == 3) {
231           /* Tetrahedra are inverted */
232           if (gmsh_elem[c].numNodes == 4) {
233             PetscInt tmp = pcone[0];
234             pcone[0] = pcone[1];
235             pcone[1] = tmp;
236           }
237           /* Hexahedra are inverted */
238           if (gmsh_elem[c].numNodes == 8) {
239             PetscInt tmp = pcone[1];
240             pcone[1] = pcone[3];
241             pcone[3] = tmp;
242           }
243         }
244         ierr = DMPlexSetCone(*dm, cell, pcone);CHKERRQ(ierr);
245         cell++;
246       }
247     }
248   }
249   ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr);
250   ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
251   ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
252   ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
253   if (interpolate) {
254     DM idm = NULL;
255 
256     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
257     ierr = DMDestroy(dm);CHKERRQ(ierr);
258     *dm  = idm;
259   }
260 
261   if (!rank) {
262     /* Apply boundary IDs by finding the relevant facets with vertex joins */
263     PetscInt pcone[8], corner, vStart, vEnd;
264 
265     ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
266     for (c = 0; c < numCells; ++c) {
267       if (gmsh_elem[c].dim == dim-1) {
268         PetscInt joinSize;
269         const PetscInt *join;
270         for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) {
271           pcone[corner] = gmsh_elem[c].nodes[corner] + vStart - 1;
272         }
273         ierr = DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr);
274         if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id);
275         ierr = DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);CHKERRQ(ierr);
276         ierr = DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr);
277       }
278     }
279 
280     /* Create cell sets */
281     for (cell = 0, c = 0; c < numCells; ++c) {
282       if (gmsh_elem[c].dim == dim) {
283         if (gmsh_elem[c].numTags > 0) {
284           ierr = DMSetLabelValue(*dm, "Cell Sets", cell, gmsh_elem[c].tags[0]);CHKERRQ(ierr);
285           cell++;
286         }
287       }
288     }
289 
290     /* Create vertex sets */
291     for (c = 0; c < numCells; ++c) {
292       if (gmsh_elem[c].dim == 0) {
293         if (gmsh_elem[c].numTags > 0) {
294           ierr = DMSetLabelValue(*dm, "Vertex Sets", gmsh_elem[c].nodes[0] + vStart - 1, gmsh_elem[c].tags[0]);CHKERRQ(ierr);
295         }
296       }
297     }
298   }
299 
300   /* Read coordinates */
301   ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
302   ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
303   ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr);
304   ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr);
305   for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
306     ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr);
307     ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr);
308   }
309   ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
310   ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
311   ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
312   ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
313   ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
314   ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr);
315   ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr);
316   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
317   if (!rank) {
318     for (v = 0; v < numVertices; ++v) {
319       for (d = 0; d < dim; ++d) {
320         coords[v*dim+d] = coordsIn[v*3+d];
321       }
322     }
323   }
324   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
325   ierr = PetscFree(coordsIn);CHKERRQ(ierr);
326   ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr);
327   ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
328   /* Clean up intermediate storage */
329   if (!rank || binary) ierr = PetscFree(gmsh_elem);CHKERRQ(ierr);
330   ierr = PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr);
331   PetscFunctionReturn(0);
332 }
333 
334 #undef __FUNCT__
335 #define __FUNCT__ "DMPlexCreateGmsh_ReadElement"
336 PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscInt numCells, PetscBool binary, PetscBool byteSwap, GmshElement **gmsh_elems)
337 {
338   PetscInt       c, p;
339   GmshElement   *elements;
340   int            i, cellType, dim, numNodes, numElem, numTags;
341   int            ibuf[16];
342   PetscErrorCode ierr;
343 
344   PetscFunctionBegin;
345   ierr = PetscMalloc1(numCells, &elements);CHKERRQ(ierr);
346   for (c = 0; c < numCells;) {
347     ierr = PetscViewerRead(viewer, &ibuf, 3, NULL, PETSC_ENUM);CHKERRQ(ierr);
348     if (byteSwap) ierr = PetscByteSwap(&ibuf, PETSC_ENUM, 3);CHKERRQ(ierr);
349     if (binary) {
350       cellType = ibuf[0];
351       numElem = ibuf[1];
352       numTags = ibuf[2];
353     } else {
354       elements[c].id = ibuf[0];
355       cellType = ibuf[1];
356       numTags = ibuf[2];
357       numElem = 1;
358     }
359     switch (cellType) {
360     case 1: /* 2-node line */
361       dim = 1;
362       numNodes = 2;
363       break;
364     case 2: /* 3-node triangle */
365       dim = 2;
366       numNodes = 3;
367       break;
368     case 3: /* 4-node quadrangle */
369       dim = 2;
370       numNodes = 4;
371       break;
372     case 4: /* 4-node tetrahedron */
373       dim  = 3;
374       numNodes = 4;
375       break;
376     case 5: /* 8-node hexahedron */
377       dim = 3;
378       numNodes = 8;
379       break;
380     case 15: /* 1-node vertex */
381       dim = 0;
382       numNodes = 1;
383       break;
384     default:
385       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
386     }
387     if (binary) {
388       const PetscInt nint = numNodes + numTags + 1;
389       for (i = 0; i < numElem; ++i, ++c) {
390         /* Loop over inner binary element block */
391         elements[c].dim = dim;
392         elements[c].numNodes = numNodes;
393         elements[c].numTags = numTags;
394 
395         ierr = PetscViewerRead(viewer, &ibuf, nint, NULL, PETSC_ENUM);CHKERRQ(ierr);
396         if (byteSwap) ierr = PetscByteSwap( &ibuf, PETSC_ENUM, nint);CHKERRQ(ierr);
397         elements[c].id = ibuf[0];
398         for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p];
399         for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
400       }
401     } else {
402       elements[c].dim = dim;
403       elements[c].numNodes = numNodes;
404       elements[c].numTags = numTags;
405       ierr = PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);CHKERRQ(ierr);
406       ierr = PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);CHKERRQ(ierr);
407       c++;
408     }
409   }
410   *gmsh_elems = elements;
411   PetscFunctionReturn(0);
412 }
413