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