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