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