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