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; 190 PetscReal *coordsIn = NULL; 191 PetscInt dim = 0, coordSize, c, v, d, r, cell; 192 int i, numVertices = 0, numCells = 0, trueNumCells = 0, numRegions = 0, snum; 193 PetscMPIInt num_proc, rank; 194 char line[PETSC_MAX_PATH_LEN]; 195 PetscBool match, binary, bswap = PETSC_FALSE; 196 PetscErrorCode ierr; 197 198 PetscFunctionBegin; 199 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 200 ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr); 201 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 202 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 203 ierr = PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 204 ierr = PetscViewerGetType(viewer, &vtype);CHKERRQ(ierr); 205 ierr = PetscStrcmp(vtype, PETSCVIEWERBINARY, &binary);CHKERRQ(ierr); 206 if (!rank || binary) { 207 PetscBool match; 208 int fileType, dataSize; 209 float version; 210 211 /* Read header */ 212 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 213 ierr = PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 214 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 215 ierr = PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);CHKERRQ(ierr); 216 snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize); 217 if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line); 218 if (version < 2.0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file must be at least version 2.0"); 219 if (dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize); 220 if (binary) { 221 int checkInt; 222 ierr = PetscViewerRead(viewer, &checkInt, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 223 if (checkInt != 1) { 224 ierr = PetscByteSwap(&checkInt, PETSC_ENUM, 1);CHKERRQ(ierr); 225 if (checkInt == 1) bswap = PETSC_TRUE; 226 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh binary file", fileType); 227 } 228 } else if (fileType) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File type %d is not a valid Gmsh ASCII file", fileType); 229 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 230 ierr = PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 231 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 232 /* OPTIONAL Read physical names */ 233 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 234 ierr = PetscStrncmp(line, "$PhysicalNames", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 235 if (match) { 236 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 237 snum = sscanf(line, "%d", &numRegions); 238 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 239 for (r = 0; r < numRegions; ++r) {ierr = PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);CHKERRQ(ierr);} 240 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 241 ierr = PetscStrncmp(line, "$EndPhysicalNames", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 242 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 243 /* Initial read for vertex section */ 244 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 245 } 246 /* Read vertices */ 247 ierr = PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 248 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 249 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 250 snum = sscanf(line, "%d", &numVertices); 251 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 252 ierr = PetscMalloc1(numVertices*3, &coordsIn);CHKERRQ(ierr); 253 if (binary) { 254 size_t doubleSize, intSize; 255 PetscInt elementSize; 256 char *buffer; 257 double *baseptr; 258 ierr = PetscDataTypeGetSize(PETSC_ENUM, &intSize);CHKERRQ(ierr); 259 ierr = PetscDataTypeGetSize(PETSC_DOUBLE, &doubleSize);CHKERRQ(ierr); 260 elementSize = (intSize + 3*doubleSize); 261 ierr = PetscMalloc1(elementSize*numVertices, &buffer);CHKERRQ(ierr); 262 ierr = PetscViewerRead(viewer, buffer, elementSize*numVertices, NULL, PETSC_CHAR);CHKERRQ(ierr); 263 if (bswap) ierr = PetscByteSwap(buffer, PETSC_CHAR, elementSize*numVertices);CHKERRQ(ierr); 264 for (v = 0; v < numVertices; ++v) { 265 baseptr = ((double*)(buffer+v*elementSize+intSize)); 266 coordsIn[v*3+0] = (PetscReal) baseptr[0]; 267 coordsIn[v*3+1] = (PetscReal) baseptr[1]; 268 coordsIn[v*3+2] = (PetscReal) baseptr[2]; 269 } 270 ierr = PetscFree(buffer);CHKERRQ(ierr); 271 } else { 272 for (v = 0; v < numVertices; ++v) { 273 ierr = PetscViewerRead(viewer, &i, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 274 ierr = PetscViewerRead(viewer, &(coordsIn[v*3]), 3, NULL, PETSC_REAL);CHKERRQ(ierr); 275 if (i != (int)v+1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid node number %D should be %D", i, v+1); 276 } 277 } 278 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 279 ierr = PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 280 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 281 /* Read cells */ 282 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 283 ierr = PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 284 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 285 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 286 snum = sscanf(line, "%d", &numCells); 287 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 288 } 289 290 if (!rank || binary) { 291 /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the 292 file contents multiple times to figure out the true number of cells and facets 293 in the given mesh. To make this more efficient we read the file contents only 294 once and store them in memory, while determining the true number of cells. */ 295 ierr = DMPlexCreateGmsh_ReadElement(viewer, numCells, binary, bswap, &gmsh_elem);CHKERRQ(ierr); 296 for (trueNumCells=0, c = 0; c < numCells; ++c) { 297 if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;} 298 if (gmsh_elem[c].dim == dim) trueNumCells++; 299 } 300 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 301 ierr = PetscStrncmp(line, "$EndElements", 12, &match);CHKERRQ(ierr); 302 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 303 } 304 /* For binary we read on all ranks, but only build the plex on rank 0 */ 305 if (binary && rank) {trueNumCells = 0; numVertices = 0;}; 306 /* Allocate the cell-vertex mesh */ 307 ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr); 308 if (!rank) { 309 for (cell = 0, c = 0; c < numCells; ++c) { 310 if (gmsh_elem[c].dim == dim) { 311 ierr = DMPlexSetConeSize(*dm, cell, gmsh_elem[c].numNodes);CHKERRQ(ierr); 312 cell++; 313 } 314 } 315 } 316 ierr = DMSetUp(*dm);CHKERRQ(ierr); 317 /* Add cell-vertex connections */ 318 if (!rank) { 319 PetscInt pcone[8], corner; 320 for (cell = 0, c = 0; c < numCells; ++c) { 321 if (gmsh_elem[c].dim == dim) { 322 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 323 pcone[corner] = gmsh_elem[c].nodes[corner] + trueNumCells-1; 324 } 325 if (dim == 3) { 326 /* Tetrahedra are inverted */ 327 if (gmsh_elem[c].numNodes == 4) { 328 PetscInt tmp = pcone[0]; 329 pcone[0] = pcone[1]; 330 pcone[1] = tmp; 331 } 332 /* Hexahedra are inverted */ 333 if (gmsh_elem[c].numNodes == 8) { 334 PetscInt tmp = pcone[1]; 335 pcone[1] = pcone[3]; 336 pcone[3] = tmp; 337 } 338 } 339 if (dim == 2) { 340 /* Triangles are inverted */ 341 if (gmsh_elem[c].numNodes == 3) { 342 PetscInt tmp = pcone[0]; 343 pcone[0] = pcone[1]; 344 pcone[1] = tmp; 345 } 346 } 347 ierr = DMPlexSetCone(*dm, cell, pcone);CHKERRQ(ierr); 348 cell++; 349 } 350 } 351 } 352 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 353 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 354 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 355 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 356 if (interpolate) { 357 DM idm = NULL; 358 359 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 360 ierr = DMDestroy(dm);CHKERRQ(ierr); 361 *dm = idm; 362 } 363 364 if (!rank) { 365 /* Apply boundary IDs by finding the relevant facets with vertex joins */ 366 PetscInt pcone[8], corner, vStart, vEnd; 367 368 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 369 for (c = 0; c < numCells; ++c) { 370 if (gmsh_elem[c].dim == dim-1) { 371 PetscInt joinSize; 372 const PetscInt *join; 373 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 374 pcone[corner] = gmsh_elem[c].nodes[corner] + vStart - 1; 375 } 376 ierr = DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 377 if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for element %d", gmsh_elem[c].id); 378 ierr = DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);CHKERRQ(ierr); 379 ierr = DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 380 } 381 } 382 383 /* Create cell sets */ 384 for (cell = 0, c = 0; c < numCells; ++c) { 385 if (gmsh_elem[c].dim == dim) { 386 if (gmsh_elem[c].numTags > 0) { 387 ierr = DMSetLabelValue(*dm, "Cell Sets", cell, gmsh_elem[c].tags[0]);CHKERRQ(ierr); 388 cell++; 389 } 390 } 391 } 392 393 /* Create vertex sets */ 394 for (c = 0; c < numCells; ++c) { 395 if (gmsh_elem[c].dim == 0) { 396 if (gmsh_elem[c].numTags > 0) { 397 ierr = DMSetLabelValue(*dm, "Vertex Sets", gmsh_elem[c].nodes[0] + vStart - 1, gmsh_elem[c].tags[0]);CHKERRQ(ierr); 398 } 399 } 400 } 401 } 402 403 /* Read coordinates */ 404 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 405 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 406 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 407 ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr); 408 for (v = trueNumCells; v < trueNumCells+numVertices; ++v) { 409 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 410 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 411 } 412 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 413 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 414 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 415 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 416 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 417 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr); 418 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 419 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 420 if (!rank) { 421 for (v = 0; v < numVertices; ++v) { 422 for (d = 0; d < dim; ++d) { 423 coords[v*dim+d] = coordsIn[v*3+d]; 424 } 425 } 426 } 427 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 428 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 429 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 430 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 431 /* Clean up intermediate storage */ 432 if (!rank || binary) ierr = PetscFree(gmsh_elem);CHKERRQ(ierr); 433 ierr = PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 434 PetscFunctionReturn(0); 435 } 436