xref: /petsc/src/dm/impls/plex/plexgmsh.c (revision 68f0caa67e2a15036f113eae346582c0146a1cbf)
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 
297   /* Read coordinates */
298   ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
299   ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
300   ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr);
301   ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr);
302   for (v = trueNumCells; v < trueNumCells+numVertices; ++v) {
303     ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr);
304     ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr);
305   }
306   ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
307   ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
308   ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
309   ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
310   ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
311   ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr);
312   ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr);
313   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
314   if (!rank) {
315     for (v = 0; v < numVertices; ++v) {
316       for (d = 0; d < dim; ++d) {
317         coords[v*dim+d] = coordsIn[v*3+d];
318       }
319     }
320   }
321   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
322   ierr = PetscFree(coordsIn);CHKERRQ(ierr);
323   ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr);
324   ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
325   /* Clean up intermediate storage */
326   if (!rank || binary) ierr = PetscFree(gmsh_elem);CHKERRQ(ierr);
327   ierr = PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr);
328   PetscFunctionReturn(0);
329 }
330 
331 #undef __FUNCT__
332 #define __FUNCT__ "DMPlexCreateGmsh_ReadElement"
333 PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscInt numCells, PetscBool binary, PetscBool byteSwap, GmshElement **gmsh_elems)
334 {
335   PetscInt       c, p;
336   GmshElement   *elements;
337   int            i, cellType, dim, numNodes, numElem, numTags;
338   int            ibuf[16];
339   PetscErrorCode ierr;
340 
341   PetscFunctionBegin;
342   ierr = PetscMalloc1(numCells, &elements);CHKERRQ(ierr);
343   for (c = 0; c < numCells;) {
344     ierr = PetscViewerRead(viewer, &ibuf, 3, NULL, PETSC_ENUM);CHKERRQ(ierr);
345     if (byteSwap) ierr = PetscByteSwap(&ibuf, PETSC_ENUM, 3);CHKERRQ(ierr);
346     if (binary) {
347       cellType = ibuf[0];
348       numElem = ibuf[1];
349       numTags = ibuf[2];
350     } else {
351       elements[c].id = ibuf[0];
352       cellType = ibuf[1];
353       numTags = ibuf[2];
354       numElem = 1;
355     }
356     switch (cellType) {
357     case 1: /* 2-node line */
358       dim = 1;
359       numNodes = 2;
360       break;
361     case 2: /* 3-node triangle */
362       dim = 2;
363       numNodes = 3;
364       break;
365     case 3: /* 4-node quadrangle */
366       dim = 2;
367       numNodes = 4;
368       break;
369     case 4: /* 4-node tetrahedron */
370       dim  = 3;
371       numNodes = 4;
372       break;
373     case 5: /* 8-node hexahedron */
374       dim = 3;
375       numNodes = 8;
376       break;
377     case 15: /* 1-node vertex */
378       dim = 0;
379       numNodes = 1;
380       break;
381     default:
382       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType);
383     }
384     if (binary) {
385       const PetscInt nint = numNodes + numTags + 1;
386       for (i = 0; i < numElem; ++i, ++c) {
387         /* Loop over inner binary element block */
388         elements[c].dim = dim;
389         elements[c].numNodes = numNodes;
390         elements[c].numTags = numTags;
391 
392         ierr = PetscViewerRead(viewer, &ibuf, nint, NULL, PETSC_ENUM);CHKERRQ(ierr);
393         if (byteSwap) ierr = PetscByteSwap( &ibuf, PETSC_ENUM, nint);CHKERRQ(ierr);
394         elements[c].id = ibuf[0];
395         for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p];
396         for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p];
397       }
398     } else {
399       elements[c].dim = dim;
400       elements[c].numNodes = numNodes;
401       elements[c].numTags = numTags;
402       ierr = PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);CHKERRQ(ierr);
403       ierr = PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);CHKERRQ(ierr);
404       c++;
405     }
406   }
407   *gmsh_elems = elements;
408   PetscFunctionReturn(0);
409 }
410