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; 92 PetscReal *coordsIn = NULL; 93 PetscInt dim = 0, coordSize, c, v, d, r, cell; 94 int i, numVertices = 0, numCells = 0, trueNumCells = 0, numRegions = 0, snum; 95 PetscMPIInt num_proc, rank; 96 char line[PETSC_MAX_PATH_LEN]; 97 PetscBool match, binary, bswap = PETSC_FALSE; 98 PetscErrorCode ierr; 99 100 PetscFunctionBegin; 101 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 102 ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr); 103 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 104 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 105 ierr = PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 106 ierr = PetscViewerGetType(viewer, &vtype);CHKERRQ(ierr); 107 ierr = PetscStrcmp(vtype, PETSCVIEWERBINARY, &binary);CHKERRQ(ierr); 108 if (!rank || binary) { 109 PetscBool match; 110 int fileType, dataSize; 111 float version; 112 113 /* Read header */ 114 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 115 ierr = PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 116 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 117 ierr = PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);CHKERRQ(ierr); 118 snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize); 119 if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line); 120 if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0"); 121 if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize); 122 if (binary) { 123 int checkInt; 124 ierr = PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 125 if (checkInt != 1) { 126 ierr = PetscByteSwap(&checkInt, PETSC_ENUM, 1);CHKERRQ(ierr); 127 if (checkInt == 1) bswap = PETSC_TRUE; 128 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType); 129 } 130 } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType); 131 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 132 ierr = PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 133 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 134 /* OPTIONAL Read physical names */ 135 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 136 ierr = PetscStrncmp(line, "$PhysicalNames", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 137 if (match) { 138 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 139 snum = sscanf(line, "%d", &numRegions); 140 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 141 for (r = 0; r < numRegions; ++r) {ierr = PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);CHKERRQ(ierr);} 142 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 143 ierr = PetscStrncmp(line, "$EndPhysicalNames", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 144 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 145 /* Initial read for vertex section */ 146 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 147 } 148 /* Read vertices */ 149 ierr = PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 150 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 151 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 152 snum = sscanf(line, "%d", &numVertices); 153 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 154 ierr = PetscMalloc1(numVertices*3, &coordsIn);CHKERRQ(ierr); 155 if (binary) { 156 size_t doubleSize, intSize; 157 PetscInt elementSize; 158 char *buffer; 159 PetscReal *baseptr; 160 ierr = PetscDataTypeGetSize(PETSC_ENUM, &intSize);CHKERRQ(ierr); 161 ierr = PetscDataTypeGetSize(PETSC_DOUBLE, &doubleSize);CHKERRQ(ierr); 162 elementSize = (intSize + 3*doubleSize); 163 ierr = PetscMalloc1(elementSize*numVertices, &buffer);CHKERRQ(ierr); 164 ierr = PetscViewerRead(viewer, buffer, elementSize*numVertices, NULL, PETSC_CHAR);CHKERRQ(ierr); 165 if (bswap) ierr = PetscByteSwap(buffer, PETSC_CHAR, elementSize*numVertices);CHKERRQ(ierr); 166 for (v = 0; v < numVertices; ++v) { 167 baseptr = ((PetscReal*)(buffer+v*elementSize+intSize)); 168 coordsIn[v*3+0] = baseptr[0]; 169 coordsIn[v*3+1] = baseptr[1]; 170 coordsIn[v*3+2] = baseptr[2]; 171 } 172 ierr = PetscFree(buffer);CHKERRQ(ierr); 173 } else { 174 for (v = 0; v < numVertices; ++v) { 175 ierr = PetscViewerRead(viewer, &i, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 176 ierr = PetscViewerRead(viewer, &(coordsIn[v*3]), 3, NULL, PETSC_REAL);CHKERRQ(ierr); 177 if (i != (int)v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %D should be %D", i, v+1); 178 } 179 } 180 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 181 ierr = PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 182 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 183 /* Read cells */ 184 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 185 ierr = PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 186 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 187 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 188 snum = sscanf(line, "%d", &numCells); 189 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 190 } 191 192 if (!rank || binary) { 193 /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the 194 file contents multiple times to figure out the true number of cells and facets 195 in the given mesh. To make this more efficient we read the file contents only 196 once and store them in memory, while determining the true number of cells. */ 197 ierr = DMPlexCreateGmsh_ReadElement(viewer, numCells, binary, bswap, &gmsh_elem);CHKERRQ(ierr); 198 for (trueNumCells=0, c = 0; c < numCells; ++c) { 199 if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;} 200 if (gmsh_elem[c].dim == dim) trueNumCells++; 201 } 202 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 203 ierr = PetscStrncmp(line, "$EndElements", 12, &match);CHKERRQ(ierr); 204 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 205 } 206 /* For binary we read on all ranks, but only build the plex on rank 0 */ 207 if (binary && rank) {trueNumCells = 0; numVertices = 0;}; 208 /* Allocate the cell-vertex mesh */ 209 ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr); 210 if (!rank) { 211 for (cell = 0, c = 0; c < numCells; ++c) { 212 if (gmsh_elem[c].dim == dim) { 213 ierr = DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);CHKERRQ(ierr); 214 cell++; 215 } 216 } 217 } 218 ierr = DMSetUp(*dm);CHKERRQ(ierr); 219 /* Add cell-vertex connections */ 220 if (!rank) { 221 PetscInt pcone[8], corner; 222 for (cell = 0, c = 0; c < numCells; ++c) { 223 if (gmsh_elem[c].dim == dim) { 224 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 225 pcone[corner] = gmsh_elem[c].nodes[corner] + trueNumCells-1; 226 } 227 if (dim == 3) { 228 /* Tetrahedra are inverted */ 229 if (gmsh_elem[c].numNodes == 4) { 230 PetscInt tmp = pcone[0]; 231 pcone[0] = pcone[1]; 232 pcone[1] = tmp; 233 } 234 /* Hexahedra are inverted */ 235 if (gmsh_elem[c].numNodes == 8) { 236 PetscInt tmp = pcone[1]; 237 pcone[1] = pcone[3]; 238 pcone[3] = tmp; 239 } 240 } 241 ierr = DMPlexSetCone(*dm, cell, pcone);CHKERRQ(ierr); 242 cell++; 243 } 244 } 245 } 246 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 247 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 248 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 249 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 250 if (interpolate) { 251 DM idm = NULL; 252 253 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 254 ierr = DMDestroy(dm);CHKERRQ(ierr); 255 *dm = idm; 256 } 257 258 if (!rank) { 259 /* Apply boundary IDs by finding the relevant facets with vertex joins */ 260 PetscInt pcone[8], corner, vStart, vEnd; 261 262 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 263 for (c = 0; c < numCells; ++c) { 264 if (gmsh_elem[c].dim == dim-1) { 265 PetscInt joinSize; 266 const PetscInt *join; 267 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 268 pcone[corner] = gmsh_elem[c].nodes[corner] + vStart - 1; 269 } 270 ierr = DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 271 if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id); 272 ierr = DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);CHKERRQ(ierr); 273 ierr = DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 274 } 275 } 276 277 /* Create cell sets */ 278 for (cell = 0, c = 0; c < numCells; ++c) { 279 if (gmsh_elem[c].dim == dim) { 280 if (gmsh_elem[c].numTags > 0) { 281 ierr = DMSetLabelValue(*dm, "Cell Sets", cell, gmsh_elem[c].tags[0]);CHKERRQ(ierr); 282 cell++; 283 } 284 } 285 } 286 287 /* Create vertex sets */ 288 for (c = 0; c < numCells; ++c) { 289 if (gmsh_elem[c].dim == 0) { 290 if (gmsh_elem[c].numTags > 0) { 291 ierr = DMSetLabelValue(*dm, "Vertex Sets", gmsh_elem[c].nodes[0] + vStart - 1, gmsh_elem[c].tags[0]);CHKERRQ(ierr); 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 PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer viewer, PetscInt numCells, PetscBool binary, PetscBool byteSwap, GmshElement **gmsh_elems) 332 { 333 PetscInt c, p; 334 GmshElement *elements; 335 int i, cellType, dim, numNodes, numElem, numTags; 336 int ibuf[16]; 337 PetscErrorCode ierr; 338 339 PetscFunctionBegin; 340 ierr = PetscMalloc1(numCells, &elements);CHKERRQ(ierr); 341 for (c = 0; c < numCells;) { 342 ierr = PetscViewerRead(viewer, &ibuf, 3, NULL, PETSC_ENUM);CHKERRQ(ierr); 343 if (byteSwap) ierr = PetscByteSwap(&ibuf, PETSC_ENUM, 3);CHKERRQ(ierr); 344 if (binary) { 345 cellType = ibuf[0]; 346 numElem = ibuf[1]; 347 numTags = ibuf[2]; 348 } else { 349 elements[c].id = ibuf[0]; 350 cellType = ibuf[1]; 351 numTags = ibuf[2]; 352 numElem = 1; 353 } 354 switch (cellType) { 355 case 1: /* 2-node line */ 356 dim = 1; 357 numNodes = 2; 358 break; 359 case 2: /* 3-node triangle */ 360 dim = 2; 361 numNodes = 3; 362 break; 363 case 3: /* 4-node quadrangle */ 364 dim = 2; 365 numNodes = 4; 366 break; 367 case 4: /* 4-node tetrahedron */ 368 dim = 3; 369 numNodes = 4; 370 break; 371 case 5: /* 8-node hexahedron */ 372 dim = 3; 373 numNodes = 8; 374 break; 375 case 15: /* 1-node vertex */ 376 dim = 0; 377 numNodes = 1; 378 break; 379 default: 380 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 381 } 382 if (binary) { 383 const PetscInt nint = numNodes + numTags + 1; 384 for (i = 0; i < numElem; ++i, ++c) { 385 /* Loop over inner binary element block */ 386 elements[c].dim = dim; 387 elements[c].numNodes = numNodes; 388 elements[c].numTags = numTags; 389 390 ierr = PetscViewerRead(viewer, &ibuf, nint, NULL, PETSC_ENUM);CHKERRQ(ierr); 391 if (byteSwap) ierr = PetscByteSwap( &ibuf, PETSC_ENUM, nint);CHKERRQ(ierr); 392 elements[c].id = ibuf[0]; 393 for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p]; 394 for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p]; 395 } 396 } else { 397 elements[c].dim = dim; 398 elements[c].numNodes = numNodes; 399 elements[c].numTags = numTags; 400 ierr = PetscViewerRead(viewer, elements[c].tags, elements[c].numTags, NULL, PETSC_ENUM);CHKERRQ(ierr); 401 ierr = PetscViewerRead(viewer, elements[c].nodes, elements[c].numNodes, NULL, PETSC_ENUM);CHKERRQ(ierr); 402 c++; 403 } 404 } 405 *gmsh_elems = elements; 406 PetscFunctionReturn(0); 407 } 408