1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petsc/private/hashmapi.h> 4 5 /*@C 6 DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file 7 8 + comm - The MPI communicator 9 . filename - Name of the Gmsh file 10 - interpolate - Create faces and edges in the mesh 11 12 Output Parameter: 13 . dm - The DM object representing the mesh 14 15 Level: beginner 16 17 .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 18 @*/ 19 PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 20 { 21 PetscViewer viewer; 22 PetscMPIInt rank; 23 int fileType; 24 PetscViewerType vtype; 25 PetscErrorCode ierr; 26 27 PetscFunctionBegin; 28 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 29 30 /* Determine Gmsh file type (ASCII or binary) from file header */ 31 if (!rank) { 32 PetscViewer vheader; 33 char line[PETSC_MAX_PATH_LEN]; 34 PetscBool match; 35 int snum; 36 float version; 37 38 ierr = PetscViewerCreate(PETSC_COMM_SELF, &vheader);CHKERRQ(ierr); 39 ierr = PetscViewerSetType(vheader, PETSCVIEWERASCII);CHKERRQ(ierr); 40 ierr = PetscViewerFileSetMode(vheader, FILE_MODE_READ);CHKERRQ(ierr); 41 ierr = PetscViewerFileSetName(vheader, filename);CHKERRQ(ierr); 42 /* Read only the first two lines of the Gmsh file */ 43 ierr = PetscViewerRead(vheader, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 44 ierr = PetscStrncmp(line, "$MeshFormat", sizeof(line), &match);CHKERRQ(ierr); 45 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 46 ierr = PetscViewerRead(vheader, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 47 snum = sscanf(line, "%f %d", &version, &fileType); 48 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line); 49 if (version < 2.2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at least 2.2", (double)version); 50 if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f not supported", (double)version); 51 if (version > 4.1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at most 4.1", (double)version); 52 ierr = PetscViewerDestroy(&vheader);CHKERRQ(ierr); 53 } 54 ierr = MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);CHKERRQ(ierr); 55 vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY; 56 57 /* Create appropriate viewer and build plex */ 58 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 59 ierr = PetscViewerSetType(viewer, vtype);CHKERRQ(ierr); 60 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 61 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 62 ierr = DMPlexCreateGmsh(comm, viewer, interpolate, dm);CHKERRQ(ierr); 63 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 64 PetscFunctionReturn(0); 65 } 66 67 typedef struct { 68 PetscViewer viewer; 69 int fileFormat; 70 int dataSize; 71 PetscBool binary; 72 PetscBool byteSwap; 73 size_t wlen; 74 void *wbuf; 75 size_t slen; 76 void *sbuf; 77 } GmshFile; 78 79 static PetscErrorCode GmshBufferGet(GmshFile *gmsh, size_t count, size_t eltsize, void *buf) 80 { 81 size_t size = count * eltsize; 82 PetscErrorCode ierr; 83 84 PetscFunctionBegin; 85 if (gmsh->wlen < size) { 86 ierr = PetscFree(gmsh->wbuf);CHKERRQ(ierr); 87 ierr = PetscMalloc(size, &gmsh->wbuf);CHKERRQ(ierr); 88 gmsh->wlen = size; 89 } 90 *(void**)buf = size ? gmsh->wbuf : NULL; 91 PetscFunctionReturn(0); 92 } 93 94 static PetscErrorCode GmshBufferSizeGet(GmshFile *gmsh, size_t count, void *buf) 95 { 96 size_t dataSize = (size_t)gmsh->dataSize; 97 size_t size = count * dataSize; 98 PetscErrorCode ierr; 99 100 PetscFunctionBegin; 101 if (gmsh->slen < size) { 102 ierr = PetscFree(gmsh->sbuf);CHKERRQ(ierr); 103 ierr = PetscMalloc(size, &gmsh->sbuf);CHKERRQ(ierr); 104 gmsh->slen = size; 105 } 106 *(void**)buf = size ? gmsh->sbuf : NULL; 107 PetscFunctionReturn(0); 108 } 109 110 static PetscErrorCode GmshRead(GmshFile *gmsh, void *buf, PetscInt count, PetscDataType dtype) 111 { 112 PetscErrorCode ierr; 113 PetscFunctionBegin; 114 ierr = PetscViewerRead(gmsh->viewer, buf, count, NULL, dtype);CHKERRQ(ierr); 115 if (gmsh->byteSwap) {ierr = PetscByteSwap(buf, dtype, count);CHKERRQ(ierr);} 116 PetscFunctionReturn(0); 117 } 118 119 static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscInt count) 120 { 121 PetscErrorCode ierr; 122 PetscFunctionBegin; 123 ierr = PetscViewerRead(gmsh->viewer, buf, count, NULL, PETSC_STRING);CHKERRQ(ierr); 124 PetscFunctionReturn(0); 125 } 126 127 static PetscErrorCode GmshReadSize(GmshFile *gmsh, PetscInt *buf, PetscInt count) 128 { 129 PetscInt i; 130 size_t dataSize = (size_t)gmsh->dataSize; 131 PetscErrorCode ierr; 132 133 PetscFunctionBegin; 134 if (dataSize == sizeof(PetscInt)) { 135 ierr = GmshRead(gmsh, buf, count, PETSC_INT);CHKERRQ(ierr); 136 } else if (dataSize == sizeof(int)) { 137 int *ibuf = NULL; 138 ierr = GmshBufferSizeGet(gmsh, count, &ibuf); 139 ierr = GmshRead(gmsh, ibuf, count, PETSC_ENUM);CHKERRQ(ierr); 140 for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i]; 141 } else if (dataSize == sizeof(long)) { 142 long *ibuf = NULL; 143 ierr = GmshBufferSizeGet(gmsh, count, &ibuf); 144 ierr = GmshRead(gmsh, ibuf, count, PETSC_LONG);CHKERRQ(ierr); 145 for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i]; 146 } else if (dataSize == sizeof(PetscInt64)) { 147 PetscInt64 *ibuf = NULL; 148 ierr = GmshBufferSizeGet(gmsh, count, &ibuf); 149 ierr = GmshRead(gmsh, ibuf, count, PETSC_INT64);CHKERRQ(ierr); 150 for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i]; 151 } 152 PetscFunctionReturn(0); 153 } 154 155 static PetscErrorCode GmshReadInt(GmshFile *gmsh, int *buf, PetscInt count) 156 { 157 PetscErrorCode ierr; 158 PetscFunctionBegin; 159 ierr = GmshRead(gmsh, buf, count, PETSC_ENUM);CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 static PetscErrorCode GmshReadDouble(GmshFile *gmsh, double *buf, PetscInt count) 164 { 165 PetscErrorCode ierr; 166 PetscFunctionBegin; 167 ierr = GmshRead(gmsh, buf, count, PETSC_DOUBLE);CHKERRQ(ierr); 168 PetscFunctionReturn(0); 169 } 170 171 typedef struct { 172 PetscInt id; /* Entity number */ 173 PetscInt dim; /* Entity dimension */ 174 double bbox[6]; /* Bounding box */ 175 PetscInt numTags; /* Size of tag array */ 176 int tags[4]; /* Tag array */ 177 } GmshEntity; 178 179 typedef struct { 180 GmshEntity *entity[4]; 181 PetscHMapI entityMap[4]; 182 } GmshEntities; 183 184 static PetscErrorCode GmshEntitiesCreate(PetscInt count[4], GmshEntities **entities) 185 { 186 PetscInt dim; 187 PetscErrorCode ierr; 188 189 PetscFunctionBegin; 190 ierr = PetscNew(entities);CHKERRQ(ierr); 191 for (dim = 0; dim < 4; ++dim) { 192 ierr = PetscCalloc1(count[dim], &(*entities)->entity[dim]);CHKERRQ(ierr); 193 ierr = PetscHMapICreate(&(*entities)->entityMap[dim]);CHKERRQ(ierr); 194 } 195 PetscFunctionReturn(0); 196 } 197 198 static PetscErrorCode GmshEntitiesAdd(GmshEntities *entities, PetscInt index, PetscInt dim, PetscInt eid, GmshEntity** entity) 199 { 200 PetscErrorCode ierr; 201 PetscFunctionBegin; 202 ierr = PetscHMapISet(entities->entityMap[dim], eid, index);CHKERRQ(ierr); 203 entities->entity[dim][index].dim = dim; 204 entities->entity[dim][index].id = eid; 205 if (entity) *entity = &entities->entity[dim][index]; 206 PetscFunctionReturn(0); 207 } 208 209 static PetscErrorCode GmshEntitiesGet(GmshEntities *entities, PetscInt dim, PetscInt eid, GmshEntity** entity) 210 { 211 PetscInt index; 212 PetscErrorCode ierr; 213 214 PetscFunctionBegin; 215 ierr = PetscHMapIGet(entities->entityMap[dim], eid, &index);CHKERRQ(ierr); 216 *entity = &entities->entity[dim][index]; 217 PetscFunctionReturn(0); 218 } 219 220 static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities) 221 { 222 PetscInt dim; 223 PetscErrorCode ierr; 224 225 PetscFunctionBegin; 226 if (!*entities) PetscFunctionReturn(0); 227 for (dim = 0; dim < 4; ++dim) { 228 ierr = PetscFree((*entities)->entity[dim]);CHKERRQ(ierr); 229 ierr = PetscHMapIDestroy(&(*entities)->entityMap[dim]);CHKERRQ(ierr); 230 } 231 ierr = PetscFree((*entities));CHKERRQ(ierr); 232 PetscFunctionReturn(0); 233 } 234 235 typedef struct { 236 PetscInt id; /* Entity number */ 237 PetscInt dim; /* Entity dimension */ 238 PetscInt cellType; /* Cell type */ 239 PetscInt numNodes; /* Size of node array */ 240 PetscInt nodes[8]; /* Node array */ 241 PetscInt numTags; /* Size of tag array */ 242 int tags[4]; /* Tag array */ 243 } GmshElement; 244 245 static PetscErrorCode DMPlexCreateGmsh_ReadNodes_v22(GmshFile *gmsh, int shift, PetscInt *numVertices, double **gmsh_nodes) 246 { 247 PetscViewer viewer = gmsh->viewer; 248 PetscBool byteSwap = gmsh->byteSwap; 249 char line[PETSC_MAX_PATH_LEN]; 250 int v, num, nid, snum; 251 double *coordinates; 252 PetscErrorCode ierr; 253 254 PetscFunctionBegin; 255 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 256 snum = sscanf(line, "%d", &num); 257 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 258 ierr = PetscMalloc1(num*3, &coordinates);CHKERRQ(ierr); 259 *numVertices = num; 260 *gmsh_nodes = coordinates; 261 for (v = 0; v < num; ++v) { 262 double *xyz = coordinates + v*3; 263 ierr = PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 264 if (byteSwap) {ierr = PetscByteSwap(&nid, PETSC_ENUM, 1);CHKERRQ(ierr);} 265 if (nid != v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %d", nid, v+shift); 266 ierr = PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);CHKERRQ(ierr); 267 if (byteSwap) {ierr = PetscByteSwap(xyz, PETSC_DOUBLE, 3);CHKERRQ(ierr);} 268 } 269 PetscFunctionReturn(0); 270 } 271 272 /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the 273 file contents multiple times to figure out the true number of cells and facets 274 in the given mesh. To make this more efficient we read the file contents only 275 once and store them in memory, while determining the true number of cells. */ 276 static PetscErrorCode DMPlexCreateGmsh_ReadElements_v22(GmshFile* gmsh, PETSC_UNUSED int shift, PetscInt *numCells, GmshElement **gmsh_elems) 277 { 278 PetscViewer viewer = gmsh->viewer; 279 PetscBool binary = gmsh->binary; 280 PetscBool byteSwap = gmsh->byteSwap; 281 char line[PETSC_MAX_PATH_LEN]; 282 GmshElement *elements; 283 int i, c, p, num, ibuf[1+4+512], snum; 284 int cellType, dim, numNodes, numNodesIgnore, numElem, numTags; 285 PetscErrorCode ierr; 286 287 PetscFunctionBegin; 288 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr);; 289 snum = sscanf(line, "%d", &num); 290 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 291 ierr = PetscMalloc1(num, &elements);CHKERRQ(ierr); 292 *numCells = num; 293 *gmsh_elems = elements; 294 for (c = 0; c < num;) { 295 ierr = PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);CHKERRQ(ierr); 296 if (byteSwap) {ierr = PetscByteSwap(ibuf, PETSC_ENUM, 3);CHKERRQ(ierr);} 297 if (binary) { 298 cellType = ibuf[0]; 299 numElem = ibuf[1]; 300 numTags = ibuf[2]; 301 } else { 302 elements[c].id = ibuf[0]; 303 cellType = ibuf[1]; 304 numTags = ibuf[2]; 305 numElem = 1; 306 } 307 /* http://gmsh.info/doc/texinfo/gmsh.html#MSH-ASCII-file-format */ 308 numNodesIgnore = 0; 309 switch (cellType) { 310 case 1: /* 2-node line */ 311 dim = 1; 312 numNodes = 2; 313 break; 314 case 2: /* 3-node triangle */ 315 dim = 2; 316 numNodes = 3; 317 break; 318 case 3: /* 4-node quadrangle */ 319 dim = 2; 320 numNodes = 4; 321 break; 322 case 4: /* 4-node tetrahedron */ 323 dim = 3; 324 numNodes = 4; 325 break; 326 case 5: /* 8-node hexahedron */ 327 dim = 3; 328 numNodes = 8; 329 break; 330 case 6: /* 6-node wedge */ 331 dim = 3; 332 numNodes = 6; 333 break; 334 case 8: /* 3-node 2nd order line */ 335 dim = 1; 336 numNodes = 2; 337 numNodesIgnore = 1; 338 break; 339 case 9: /* 6-node 2nd order triangle */ 340 dim = 2; 341 numNodes = 3; 342 numNodesIgnore = 3; 343 break; 344 case 13: /* 18-node 2nd wedge */ 345 dim = 3; 346 numNodes = 6; 347 numNodesIgnore = 12; 348 break; 349 case 15: /* 1-node vertex */ 350 dim = 0; 351 numNodes = 1; 352 break; 353 case 7: /* 5-node pyramid */ 354 case 10: /* 9-node 2nd order quadrangle */ 355 case 11: /* 10-node 2nd order tetrahedron */ 356 case 12: /* 27-node 2nd order hexhedron */ 357 case 14: /* 14-node 2nd order pyramid */ 358 default: 359 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 360 } 361 if (binary) { 362 const int nint = 1 + numTags + numNodes + numNodesIgnore; 363 /* Loop over element blocks */ 364 for (i = 0; i < numElem; ++i, ++c) { 365 ierr = PetscViewerRead(viewer, ibuf, nint, NULL, PETSC_ENUM);CHKERRQ(ierr); 366 if (byteSwap) {ierr = PetscByteSwap(ibuf, PETSC_ENUM, nint);CHKERRQ(ierr);} 367 elements[c].dim = dim; 368 elements[c].numNodes = numNodes; 369 elements[c].numTags = numTags; 370 elements[c].id = ibuf[0]; 371 elements[c].cellType = cellType; 372 for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[1 + p]; 373 for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[1 + numTags + p]; 374 } 375 } else { 376 const int nint = numTags + numNodes + numNodesIgnore; 377 elements[c].dim = dim; 378 elements[c].numNodes = numNodes; 379 elements[c].numTags = numTags; 380 elements[c].cellType = cellType; 381 ierr = PetscViewerRead(viewer, ibuf, nint, NULL, PETSC_ENUM);CHKERRQ(ierr); 382 for (p = 0; p < numTags; p++) elements[c].tags[p] = ibuf[p]; 383 for (p = 0; p < numNodes; p++) elements[c].nodes[p] = ibuf[numTags + p]; 384 c++; 385 } 386 } 387 PetscFunctionReturn(0); 388 } 389 390 /* 391 $Entities 392 numPoints(unsigned long) numCurves(unsigned long) 393 numSurfaces(unsigned long) numVolumes(unsigned long) 394 // points 395 tag(int) boxMinX(double) boxMinY(double) boxMinZ(double) 396 boxMaxX(double) boxMaxY(double) boxMaxZ(double) 397 numPhysicals(unsigned long) phyisicalTag[...](int) 398 ... 399 // curves 400 tag(int) boxMinX(double) boxMinY(double) boxMinZ(double) 401 boxMaxX(double) boxMaxY(double) boxMaxZ(double) 402 numPhysicals(unsigned long) physicalTag[...](int) 403 numBREPVert(unsigned long) tagBREPVert[...](int) 404 ... 405 // surfaces 406 tag(int) boxMinX(double) boxMinY(double) boxMinZ(double) 407 boxMaxX(double) boxMaxY(double) boxMaxZ(double) 408 numPhysicals(unsigned long) physicalTag[...](int) 409 numBREPCurve(unsigned long) tagBREPCurve[...](int) 410 ... 411 // volumes 412 tag(int) boxMinX(double) boxMinY(double) boxMinZ(double) 413 boxMaxX(double) boxMaxY(double) boxMaxZ(double) 414 numPhysicals(unsigned long) physicalTag[...](int) 415 numBREPSurfaces(unsigned long) tagBREPSurfaces[...](int) 416 ... 417 $EndEntities 418 */ 419 static PetscErrorCode DMPlexCreateGmsh_ReadEntities_v40(GmshFile *gmsh, GmshEntities **entities) 420 { 421 PetscViewer viewer = gmsh->viewer; 422 PetscBool byteSwap = gmsh->byteSwap; 423 long index, num, lbuf[4]; 424 int dim, eid, numTags, *ibuf, t; 425 PetscInt count[4], i; 426 GmshEntity *entity = NULL; 427 PetscErrorCode ierr; 428 429 PetscFunctionBegin; 430 ierr = PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG);CHKERRQ(ierr); 431 if (byteSwap) {ierr = PetscByteSwap(lbuf, PETSC_LONG, 4);CHKERRQ(ierr);} 432 for (i = 0; i < 4; ++i) count[i] = lbuf[i]; 433 ierr = GmshEntitiesCreate(count, entities);CHKERRQ(ierr); 434 for (dim = 0; dim < 4; ++dim) { 435 for (index = 0; index < count[dim]; ++index) { 436 ierr = PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 437 if (byteSwap) {ierr = PetscByteSwap(&eid, PETSC_ENUM, 1);CHKERRQ(ierr);} 438 ierr = GmshEntitiesAdd(*entities, (PetscInt)index, dim, eid, &entity);CHKERRQ(ierr); 439 ierr = PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE);CHKERRQ(ierr); 440 if (byteSwap) {ierr = PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6);CHKERRQ(ierr);} 441 ierr = PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 442 if (byteSwap) {ierr = PetscByteSwap(&num, PETSC_LONG, 1);CHKERRQ(ierr);} 443 ierr = GmshBufferGet(gmsh, num, sizeof(int), &ibuf);CHKERRQ(ierr); 444 ierr = PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);CHKERRQ(ierr); 445 if (byteSwap) {ierr = PetscByteSwap(ibuf, PETSC_ENUM, num);CHKERRQ(ierr);} 446 entity->numTags = numTags = (int) PetscMin(num, 4); 447 for (t = 0; t < numTags; ++t) entity->tags[t] = ibuf[t]; 448 if (dim == 0) continue; 449 ierr = PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 450 if (byteSwap) {ierr = PetscByteSwap(&num, PETSC_LONG, 1);CHKERRQ(ierr);} 451 ierr = GmshBufferGet(gmsh, num, sizeof(int), &ibuf);CHKERRQ(ierr); 452 ierr = PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);CHKERRQ(ierr); 453 if (byteSwap) {ierr = PetscByteSwap(ibuf, PETSC_ENUM, num);CHKERRQ(ierr);} 454 } 455 } 456 PetscFunctionReturn(0); 457 } 458 459 /* 460 $Nodes 461 numEntityBlocks(unsigned long) numNodes(unsigned long) 462 tagEntity(int) dimEntity(int) typeNode(int) numNodes(unsigned long) 463 tag(int) x(double) y(double) z(double) 464 ... 465 ... 466 $EndNodes 467 */ 468 static PetscErrorCode DMPlexCreateGmsh_ReadNodes_v40(GmshFile *gmsh, int shift, PetscInt *numVertices, double **gmsh_nodes) 469 { 470 PetscViewer viewer = gmsh->viewer; 471 PetscBool byteSwap = gmsh->byteSwap; 472 long block, node, v, numEntityBlocks, numTotalNodes, numNodes; 473 int info[3], nid; 474 double *coordinates; 475 char *cbuf; 476 PetscErrorCode ierr; 477 478 PetscFunctionBegin; 479 ierr = PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 480 if (byteSwap) {ierr = PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);CHKERRQ(ierr);} 481 ierr = PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 482 if (byteSwap) {ierr = PetscByteSwap(&numTotalNodes, PETSC_LONG, 1);CHKERRQ(ierr);} 483 ierr = PetscMalloc1(numTotalNodes*3, &coordinates);CHKERRQ(ierr); 484 *numVertices = numTotalNodes; 485 *gmsh_nodes = coordinates; 486 for (v = 0, block = 0; block < numEntityBlocks; ++block) { 487 ierr = PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);CHKERRQ(ierr); 488 ierr = PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 489 if (byteSwap) {ierr = PetscByteSwap(&numNodes, PETSC_LONG, 1);CHKERRQ(ierr);} 490 if (gmsh->binary) { 491 int nbytes = sizeof(int) + 3*sizeof(double); 492 ierr = GmshBufferGet(gmsh, numNodes, nbytes, &cbuf);CHKERRQ(ierr); 493 ierr = PetscViewerRead(viewer, cbuf, numNodes*nbytes, NULL, PETSC_CHAR);CHKERRQ(ierr); 494 for (node = 0; node < numNodes; ++node, ++v) { 495 char *cnid = cbuf + node*nbytes, *cxyz = cnid + sizeof(int); 496 double *xyz = coordinates + v*3; 497 #if !defined(PETSC_WORDS_BIGENDIAN) 498 ierr = PetscByteSwap(cnid, PETSC_ENUM, 1);CHKERRQ(ierr); 499 ierr = PetscByteSwap(cxyz, PETSC_DOUBLE, 3);CHKERRQ(ierr); 500 #endif 501 ierr = PetscMemcpy(&nid, cnid, sizeof(int));CHKERRQ(ierr); 502 ierr = PetscMemcpy(xyz, cxyz, 3*sizeof(double));CHKERRQ(ierr); 503 if (byteSwap) {ierr = PetscByteSwap(&nid, PETSC_ENUM, 1);CHKERRQ(ierr);} 504 if (byteSwap) {ierr = PetscByteSwap(xyz, PETSC_DOUBLE, 3);CHKERRQ(ierr);} 505 if (nid != v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %d", nid, v+shift); 506 } 507 } else { 508 for (node = 0; node < numNodes; ++node, ++v) { 509 double *xyz = coordinates + v*3; 510 ierr = PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 511 if (byteSwap) {ierr = PetscByteSwap(&nid, PETSC_ENUM, 1);CHKERRQ(ierr);} 512 if (nid != v+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %d", nid, v+shift); 513 ierr = PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);CHKERRQ(ierr); 514 if (byteSwap) {ierr = PetscByteSwap(xyz, PETSC_DOUBLE, 3);CHKERRQ(ierr);} 515 } 516 } 517 } 518 PetscFunctionReturn(0); 519 } 520 521 /* 522 $Elements 523 numEntityBlocks(unsigned long) numElements(unsigned long) 524 tagEntity(int) dimEntity(int) typeEle(int) numElements(unsigned long) 525 tag(int) numVert[...](int) 526 ... 527 ... 528 $EndElements 529 */ 530 static PetscErrorCode DMPlexCreateGmsh_ReadElements_v40(GmshFile *gmsh, PETSC_UNUSED int shift, GmshEntities *entities, PetscInt *numCells, GmshElement **gmsh_elems) 531 { 532 PetscViewer viewer = gmsh->viewer; 533 PetscBool byteSwap = gmsh->byteSwap; 534 long c, block, numEntityBlocks, numTotalElements, elem, numElements; 535 int p, info[3], *ibuf = NULL; 536 int eid, dim, numTags, *tags, cellType, numNodes; 537 GmshEntity *entity = NULL; 538 GmshElement *elements; 539 PetscErrorCode ierr; 540 541 PetscFunctionBegin; 542 ierr = PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 543 if (byteSwap) {ierr = PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);CHKERRQ(ierr);} 544 ierr = PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 545 if (byteSwap) {ierr = PetscByteSwap(&numTotalElements, PETSC_LONG, 1);CHKERRQ(ierr);} 546 ierr = PetscCalloc1(numTotalElements, &elements);CHKERRQ(ierr); 547 *numCells = numTotalElements; 548 *gmsh_elems = elements; 549 for (c = 0, block = 0; block < numEntityBlocks; ++block) { 550 ierr = PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);CHKERRQ(ierr); 551 if (byteSwap) {ierr = PetscByteSwap(info, PETSC_ENUM, 3);CHKERRQ(ierr);} 552 eid = info[0]; dim = info[1]; cellType = info[2]; 553 ierr = GmshEntitiesGet(entities, dim, eid, &entity);CHKERRQ(ierr); 554 numTags = entity->numTags; 555 tags = entity->tags; 556 switch (cellType) { 557 case 1: /* 2-node line */ 558 numNodes = 2; 559 break; 560 case 2: /* 3-node triangle */ 561 numNodes = 3; 562 break; 563 case 3: /* 4-node quadrangle */ 564 numNodes = 4; 565 break; 566 case 4: /* 4-node tetrahedron */ 567 numNodes = 4; 568 break; 569 case 5: /* 8-node hexahedron */ 570 numNodes = 8; 571 break; 572 case 6: /* 6-node wedge */ 573 numNodes = 6; 574 break; 575 case 15: /* 1-node vertex */ 576 numNodes = 1; 577 break; 578 default: 579 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 580 } 581 ierr = PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 582 if (byteSwap) {ierr = PetscByteSwap(&numElements, PETSC_LONG, 1);CHKERRQ(ierr);} 583 ierr = GmshBufferGet(gmsh, (1+numNodes)*numElements, sizeof(int), &ibuf);CHKERRQ(ierr); 584 ierr = PetscViewerRead(viewer, ibuf, (1+numNodes)*numElements, NULL, PETSC_ENUM);CHKERRQ(ierr); 585 if (byteSwap) {ierr = PetscByteSwap(ibuf, PETSC_ENUM, (1+numNodes)*numElements);CHKERRQ(ierr);} 586 for (elem = 0; elem < numElements; ++elem, ++c) { 587 int *id = ibuf + elem*(1+numNodes), *nodes = id + 1; 588 GmshElement *element = elements + c; 589 element->dim = dim; 590 element->cellType = cellType; 591 element->numNodes = numNodes; 592 element->numTags = numTags; 593 element->id = id[0]; 594 for (p = 0; p < numNodes; p++) element->nodes[p] = nodes[p]; 595 for (p = 0; p < numTags; p++) element->tags[p] = tags[p]; 596 } 597 } 598 PetscFunctionReturn(0); 599 } 600 601 static PetscErrorCode DMPlexCreateGmsh_ReadPeriodic_v40(GmshFile *gmsh, int shift, PetscInt slaveMap[], PetscBT bt) 602 { 603 PetscViewer viewer = gmsh->viewer; 604 int fileFormat = gmsh->fileFormat; 605 PetscBool binary = gmsh->binary; 606 PetscBool byteSwap = gmsh->byteSwap; 607 int numPeriodic, snum, i; 608 char line[PETSC_MAX_PATH_LEN]; 609 PetscErrorCode ierr; 610 611 PetscFunctionBegin; 612 if (fileFormat == 22 || !binary) { 613 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 614 snum = sscanf(line, "%d", &numPeriodic); 615 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 616 } else { 617 ierr = PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 618 if (byteSwap) {ierr = PetscByteSwap(&numPeriodic, PETSC_ENUM, 1);CHKERRQ(ierr);} 619 } 620 for (i = 0; i < numPeriodic; i++) { 621 int ibuf[3], slaveDim = -1, slaveTag = -1, masterTag = -1, slaveNode, masterNode; 622 long j, nNodes; 623 double affine[16]; 624 625 if (fileFormat == 22 || !binary) { 626 ierr = PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);CHKERRQ(ierr); 627 snum = sscanf(line, "%d %d %d", &slaveDim, &slaveTag, &masterTag); 628 if (snum != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 629 } else { 630 ierr = PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);CHKERRQ(ierr); 631 if (byteSwap) {ierr = PetscByteSwap(ibuf, PETSC_ENUM, 3);CHKERRQ(ierr);} 632 slaveDim = ibuf[0]; slaveTag = ibuf[1]; masterTag = ibuf[2]; 633 } 634 (void)slaveDim; (void)slaveTag; (void)masterTag; /* unused */ 635 636 if (fileFormat == 22 || !binary) { 637 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 638 snum = sscanf(line, "%ld", &nNodes); 639 if (snum != 1) { /* discard tranformation and try again */ 640 ierr = PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);CHKERRQ(ierr); 641 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 642 snum = sscanf(line, "%ld", &nNodes); 643 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 644 } 645 } else { 646 ierr = PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 647 if (byteSwap) {ierr = PetscByteSwap(&nNodes, PETSC_LONG, 1);CHKERRQ(ierr);} 648 if (nNodes == -1) { /* discard tranformation and try again */ 649 ierr = PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE);CHKERRQ(ierr); 650 ierr = PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 651 if (byteSwap) {ierr = PetscByteSwap(&nNodes, PETSC_LONG, 1);CHKERRQ(ierr);} 652 } 653 } 654 655 for (j = 0; j < nNodes; j++) { 656 if (fileFormat == 22 || !binary) { 657 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 658 snum = sscanf(line, "%d %d", &slaveNode, &masterNode); 659 if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 660 } else { 661 ierr = PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM);CHKERRQ(ierr); 662 if (byteSwap) {ierr = PetscByteSwap(ibuf, PETSC_ENUM, 2);CHKERRQ(ierr);} 663 slaveNode = ibuf[0]; masterNode = ibuf[1]; 664 } 665 slaveMap[slaveNode - shift] = masterNode - shift; 666 ierr = PetscBTSet(bt, slaveNode - shift);CHKERRQ(ierr); 667 ierr = PetscBTSet(bt, masterNode - shift);CHKERRQ(ierr); 668 } 669 } 670 PetscFunctionReturn(0); 671 } 672 673 /* 674 $Entities 675 numPoints(size_t) numCurves(size_t) 676 numSurfaces(size_t) numVolumes(size_t) 677 pointTag(int) X(double) Y(double) Z(double) 678 numPhysicalTags(size_t) physicalTag(int) ... 679 ... 680 curveTag(int) minX(double) minY(double) minZ(double) 681 maxX(double) maxY(double) maxZ(double) 682 numPhysicalTags(size_t) physicalTag(int) ... 683 numBoundingPoints(size_t) pointTag(int) ... 684 ... 685 surfaceTag(int) minX(double) minY(double) minZ(double) 686 maxX(double) maxY(double) maxZ(double) 687 numPhysicalTags(size_t) physicalTag(int) ... 688 numBoundingCurves(size_t) curveTag(int) ... 689 ... 690 volumeTag(int) minX(double) minY(double) minZ(double) 691 maxX(double) maxY(double) maxZ(double) 692 numPhysicalTags(size_t) physicalTag(int) ... 693 numBoundngSurfaces(size_t) surfaceTag(int) ... 694 ... 695 $EndEntities 696 */ 697 static PetscErrorCode DMPlexCreateGmsh_ReadEntities_v41(GmshFile *gmsh, GmshEntities **entities) 698 { 699 PetscInt count[4], index, numTags, i; 700 int dim, eid, *tags = NULL; 701 GmshEntity *entity = NULL; 702 PetscErrorCode ierr; 703 704 PetscFunctionBegin; 705 ierr = GmshReadSize(gmsh, count, 4);CHKERRQ(ierr); 706 ierr = GmshEntitiesCreate(count, entities);CHKERRQ(ierr); 707 for (dim = 0; dim < 4; ++dim) { 708 for (index = 0; index < count[dim]; ++index) { 709 ierr = GmshReadInt(gmsh, &eid, 1);CHKERRQ(ierr); 710 ierr = GmshEntitiesAdd(*entities, (PetscInt)index, dim, eid, &entity);CHKERRQ(ierr); 711 ierr = GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6);CHKERRQ(ierr); 712 ierr = GmshReadSize(gmsh, &numTags, 1);CHKERRQ(ierr); 713 ierr = GmshBufferGet(gmsh, numTags, sizeof(int), &tags);CHKERRQ(ierr); 714 ierr = GmshReadInt(gmsh, tags, numTags);CHKERRQ(ierr); 715 entity->numTags = PetscMin(numTags, 4); 716 for (i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i]; 717 if (dim == 0) continue; 718 ierr = GmshReadSize(gmsh, &numTags, 1);CHKERRQ(ierr); 719 ierr = GmshBufferGet(gmsh, numTags, sizeof(int), &tags);CHKERRQ(ierr); 720 ierr = GmshReadInt(gmsh, tags, numTags);CHKERRQ(ierr); 721 } 722 } 723 PetscFunctionReturn(0); 724 } 725 726 /* 727 $Nodes 728 numEntityBlocks(size_t) numNodes(size_t) 729 minNodeTag(size_t) maxNodeTag(size_t) 730 entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesBlock(size_t) 731 nodeTag(size_t) 732 ... 733 x(double) y(double) z(double) 734 < u(double; if parametric and entityDim = 1 or entityDim = 2) > 735 < v(double; if parametric and entityDim = 2) > 736 ... 737 ... 738 $EndNodes 739 */ 740 static PetscErrorCode DMPlexCreateGmsh_ReadNodes_v41(GmshFile *gmsh, int shift, PetscInt *numVertices, double **gmsh_nodes) 741 { 742 int info[3]; 743 PetscInt sizes[4], numEntityBlocks, numNodes, numNodesBlock = 0, *nodeTag = NULL, block, node, i; 744 double *coordinates; 745 PetscErrorCode ierr; 746 747 PetscFunctionBegin; 748 ierr = GmshReadSize(gmsh, sizes, 4);CHKERRQ(ierr); 749 numEntityBlocks = sizes[0]; numNodes = sizes[1]; 750 ierr = PetscMalloc1(numNodes*3, &coordinates);CHKERRQ(ierr); 751 *numVertices = numNodes; 752 *gmsh_nodes = coordinates; 753 for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) { 754 ierr = GmshReadInt(gmsh, info, 3);CHKERRQ(ierr); 755 ierr = GmshReadSize(gmsh, &numNodesBlock, 1);CHKERRQ(ierr); 756 ierr = GmshBufferGet(gmsh, numNodesBlock, sizeof(PetscInt), &nodeTag);CHKERRQ(ierr); 757 ierr = GmshReadSize(gmsh, nodeTag, numNodesBlock);CHKERRQ(ierr); 758 for (i = 0; i < numNodesBlock; ++i) if (nodeTag[i] != node+i+shift) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unexpected node number %d should be %d", nodeTag[i], node+i+shift); 759 if (info[2] != 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parametric coordinates not supported"); 760 ierr = GmshReadDouble(gmsh, coordinates+node*3, numNodesBlock*3);CHKERRQ(ierr); 761 } 762 PetscFunctionReturn(0); 763 } 764 765 /* 766 $Elements 767 numEntityBlocks(size_t) numElements(size_t) 768 minElementTag(size_t) maxElementTag(size_t) 769 entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t) 770 elementTag(size_t) nodeTag(size_t) ... 771 ... 772 ... 773 $EndElements 774 */ 775 static PetscErrorCode DMPlexCreateGmsh_ReadElements_v41(GmshFile *gmsh, PETSC_UNUSED int shift, GmshEntities *entities, PetscInt *numCells, GmshElement **gmsh_elems) 776 { 777 int info[3], eid, dim, cellType, *tags; 778 PetscInt sizes[4], *ibuf = NULL, numEntityBlocks, numElements, numBlockElements, numNodes, numTags, block, elem, c, p; 779 GmshEntity *entity = NULL; 780 GmshElement *elements; 781 PetscErrorCode ierr; 782 783 PetscFunctionBegin; 784 ierr = GmshReadSize(gmsh, sizes, 4);CHKERRQ(ierr); 785 numEntityBlocks = sizes[0]; numElements = sizes[1]; 786 ierr = PetscCalloc1(numElements, &elements);CHKERRQ(ierr); 787 *numCells = numElements; 788 *gmsh_elems = elements; 789 for (c = 0, block = 0; block < numEntityBlocks; ++block) { 790 ierr = GmshReadInt(gmsh, info, 3);CHKERRQ(ierr); 791 dim = info[0]; eid = info[1]; cellType = info[2]; 792 ierr = GmshEntitiesGet(entities, dim, eid, &entity);CHKERRQ(ierr); 793 numTags = entity->numTags; 794 tags = entity->tags; 795 switch (cellType) { 796 case 1: /* 2-node line */ 797 numNodes = 2; 798 break; 799 case 2: /* 3-node triangle */ 800 numNodes = 3; 801 break; 802 case 3: /* 4-node quadrangle */ 803 numNodes = 4; 804 break; 805 case 4: /* 4-node tetrahedron */ 806 numNodes = 4; 807 break; 808 case 5: /* 8-node hexahedron */ 809 numNodes = 8; 810 break; 811 case 6: /* 6-node wedge */ 812 numNodes = 6; 813 break; 814 case 15: /* 1-node vertex */ 815 numNodes = 1; 816 break; 817 default: 818 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported Gmsh element type %d", cellType); 819 } 820 ierr = GmshReadSize(gmsh, &numBlockElements, 1);CHKERRQ(ierr); 821 ierr = GmshBufferGet(gmsh, (1+numNodes)*numBlockElements, sizeof(PetscInt), &ibuf);CHKERRQ(ierr); 822 ierr = GmshReadSize(gmsh, ibuf, (1+numNodes)*numBlockElements);CHKERRQ(ierr); 823 for (elem = 0; elem < numBlockElements; ++elem, ++c) { 824 GmshElement *element = elements + c; 825 PetscInt *id = ibuf + elem*(1+numNodes), *nodes = id + 1; 826 element->id = id[0]; 827 element->dim = dim; 828 element->cellType = cellType; 829 element->numNodes = numNodes; 830 element->numTags = numTags; 831 for (p = 0; p < numNodes; p++) element->nodes[p] = nodes[p]; 832 for (p = 0; p < numTags; p++) element->tags[p] = tags[p]; 833 } 834 } 835 PetscFunctionReturn(0); 836 } 837 838 /* 839 $Periodic 840 numPeriodicLinks(size_t) 841 entityDim(int) entityTag(int) entityTagMaster(int) 842 numAffine(size_t) value(double) ... 843 numCorrespondingNodes(size_t) 844 nodeTag(size_t) nodeTagMaster(size_t) 845 ... 846 ... 847 $EndPeriodic 848 */ 849 static PetscErrorCode DMPlexCreateGmsh_ReadPeriodic_v41(GmshFile *gmsh, int shift, PetscInt slaveMap[], PetscBT bt) 850 { 851 int info[3]; 852 PetscInt numPeriodicLinks, numAffine, numCorrespondingNodes, *nodeTags = NULL, link, node; 853 double dbuf[16]; 854 PetscErrorCode ierr; 855 856 PetscFunctionBegin; 857 ierr = GmshReadSize(gmsh, &numPeriodicLinks, 1);CHKERRQ(ierr); 858 for (link = 0; link < numPeriodicLinks; ++link) { 859 ierr = GmshReadInt(gmsh, info, 3);CHKERRQ(ierr); 860 ierr = GmshReadSize(gmsh, &numAffine, 1);CHKERRQ(ierr); 861 ierr = GmshReadDouble(gmsh, dbuf, numAffine);CHKERRQ(ierr); 862 ierr = GmshReadSize(gmsh, &numCorrespondingNodes, 1);CHKERRQ(ierr); 863 ierr = GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags);CHKERRQ(ierr); 864 ierr = GmshReadSize(gmsh, nodeTags, numCorrespondingNodes*2);CHKERRQ(ierr); 865 for (node = 0; node < numCorrespondingNodes; ++node) { 866 PetscInt slaveNode = nodeTags[node*2+0] - shift; 867 PetscInt masterNode = nodeTags[node*2+1] - shift; 868 slaveMap[slaveNode] = masterNode; 869 ierr = PetscBTSet(bt, slaveNode);CHKERRQ(ierr); 870 ierr = PetscBTSet(bt, masterNode);CHKERRQ(ierr); 871 } 872 } 873 PetscFunctionReturn(0); 874 } 875 876 static PetscErrorCode DMPlexCreateGmsh_ReadMeshFormat(GmshFile *gmsh) 877 { 878 char line[PETSC_MAX_PATH_LEN]; 879 int snum, fileType, fileFormat, dataSize, checkEndian; 880 float version; 881 PetscErrorCode ierr; 882 883 PetscFunctionBegin; 884 ierr = GmshReadString(gmsh, line, 3);CHKERRQ(ierr); 885 snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize); 886 if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse Gmsh file header: %s", line); 887 if (version < 2.2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at least 2.2", (double)version); 888 if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f not supported", (double)version); 889 if (version > 4.1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Gmsh file version %3.1f must be at most 4.1", (double)version); 890 if (gmsh->binary && !fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Viewer is binary but Gmsh file is ASCII"); 891 if (!gmsh->binary && fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Viewer is ASCII but Gmsh file is binary"); 892 fileFormat = (int)(version*10); /* XXX Should use (int)roundf(version*10) ? */ 893 if (fileFormat <= 40 && dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize); 894 if (fileFormat >= 41 && dataSize != sizeof(int) && dataSize != sizeof(PetscInt64)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Data size %d is not valid for a Gmsh file", dataSize); 895 gmsh->fileFormat = fileFormat; 896 gmsh->dataSize = dataSize; 897 gmsh->byteSwap = PETSC_FALSE; 898 if (gmsh->binary) { 899 ierr = GmshReadInt(gmsh, &checkEndian, 1);CHKERRQ(ierr); 900 if (checkEndian != 1) { 901 ierr = PetscByteSwap(&checkEndian, PETSC_ENUM, 1);CHKERRQ(ierr); 902 if (checkEndian != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to detect endianness in Gmsh file header: %s", line); 903 gmsh->byteSwap = PETSC_TRUE; 904 } 905 } 906 PetscFunctionReturn(0); 907 } 908 909 static PetscErrorCode DMPlexCreateGmsh_ReadPhysicalNames(GmshFile *gmsh) 910 { 911 char line[PETSC_MAX_PATH_LEN]; 912 int snum, numRegions, region; 913 PetscErrorCode ierr; 914 915 PetscFunctionBegin; 916 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 917 snum = sscanf(line, "%d", &numRegions); 918 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 919 for (region = 0; region < numRegions; ++region) { 920 ierr = GmshReadString(gmsh, line, 3);CHKERRQ(ierr); 921 } 922 PetscFunctionReturn(0); 923 } 924 925 /*@ 926 DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer 927 928 Collective on comm 929 930 Input Parameters: 931 + comm - The MPI communicator 932 . viewer - The Viewer associated with a Gmsh file 933 - interpolate - Create faces and edges in the mesh 934 935 Output Parameter: 936 . dm - The DM object representing the mesh 937 938 Note: http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format 939 940 Level: beginner 941 942 .keywords: mesh,Gmsh 943 .seealso: DMPLEX, DMCreate() 944 @*/ 945 PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm) 946 { 947 PetscViewer parentviewer = NULL; 948 double *coordsIn = NULL; 949 GmshEntities *entities = NULL; 950 GmshElement *gmsh_elem = NULL; 951 PetscSection coordSection; 952 Vec coordinates; 953 PetscBT periodicV = NULL, periodicC = NULL; 954 PetscScalar *coords; 955 PetscInt dim = 0, embedDim = -1, coordSize, c, v, d, cell, *periodicMap = NULL, *periodicMapI = NULL, *hybridMap = NULL, cMax = PETSC_DETERMINE; 956 PetscInt numVertices = 0, numCells = 0, trueNumCells = 0; 957 int i, shift = 1; 958 PetscMPIInt rank; 959 PetscBool binary, zerobase = PETSC_FALSE, periodic = PETSC_FALSE, usemarker = PETSC_FALSE; 960 PetscBool enable_hybrid = PETSC_FALSE; 961 PetscErrorCode ierr; 962 963 PetscFunctionBegin; 964 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 965 ierr = PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_hybrid", &enable_hybrid, NULL);CHKERRQ(ierr); 966 ierr = PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_periodic", &periodic, NULL);CHKERRQ(ierr); 967 ierr = PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_use_marker", &usemarker, NULL);CHKERRQ(ierr); 968 ierr = PetscOptionsGetBool(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_zero_base", &zerobase, NULL);CHKERRQ(ierr); 969 ierr = PetscOptionsGetInt(((PetscObject) viewer)->options,((PetscObject) viewer)->prefix, "-dm_plex_gmsh_spacedim", &embedDim, NULL);CHKERRQ(ierr); 970 if (zerobase) shift = 0; 971 972 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 973 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 974 ierr = PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 975 976 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);CHKERRQ(ierr); 977 978 /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */ 979 if (binary) { 980 parentviewer = viewer; 981 ierr = PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);CHKERRQ(ierr); 982 } 983 984 if (!rank) { 985 GmshFile gmsh_, *gmsh = &gmsh_; 986 char line[PETSC_MAX_PATH_LEN]; 987 PetscBool match; 988 989 ierr = PetscMemzero(gmsh,sizeof(GmshFile));CHKERRQ(ierr); 990 gmsh->viewer = viewer; 991 gmsh->binary = binary; 992 993 /* Read mesh format */ 994 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 995 ierr = PetscStrncmp(line, "$MeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 996 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 997 ierr = DMPlexCreateGmsh_ReadMeshFormat(gmsh);CHKERRQ(ierr); 998 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 999 ierr = PetscStrncmp(line, "$EndMeshFormat", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 1000 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 1001 1002 /* OPTIONAL Read physical names */ 1003 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 1004 ierr = PetscStrncmp(line, "$PhysicalNames", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 1005 if (match) { 1006 ierr = DMPlexCreateGmsh_ReadPhysicalNames(gmsh);CHKERRQ(ierr); 1007 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 1008 ierr = PetscStrncmp(line, "$EndPhysicalNames", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 1009 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 1010 /* Initial read for entity section */ 1011 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 1012 } 1013 1014 /* Read entities */ 1015 if (gmsh->fileFormat >= 40) { 1016 ierr = PetscStrncmp(line, "$Entities", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 1017 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 1018 switch (gmsh->fileFormat) { 1019 case 41: ierr = DMPlexCreateGmsh_ReadEntities_v41(gmsh, &entities);CHKERRQ(ierr); break; 1020 default: ierr = DMPlexCreateGmsh_ReadEntities_v40(gmsh, &entities);CHKERRQ(ierr); break; 1021 } 1022 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 1023 ierr = PetscStrncmp(line, "$EndEntities", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 1024 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 1025 /* Initial read for nodes section */ 1026 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 1027 } 1028 1029 /* Read nodes */ 1030 ierr = PetscStrncmp(line, "$Nodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 1031 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 1032 switch (gmsh->fileFormat) { 1033 case 41: ierr = DMPlexCreateGmsh_ReadNodes_v41(gmsh, shift, &numVertices, &coordsIn);CHKERRQ(ierr); break; 1034 case 40: ierr = DMPlexCreateGmsh_ReadNodes_v40(gmsh, shift, &numVertices, &coordsIn);CHKERRQ(ierr); break; 1035 default: ierr = DMPlexCreateGmsh_ReadNodes_v22(gmsh, shift, &numVertices, &coordsIn);CHKERRQ(ierr); break; 1036 } 1037 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 1038 ierr = PetscStrncmp(line, "$EndNodes", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 1039 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 1040 1041 /* Read elements */ 1042 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr);; 1043 ierr = PetscStrncmp(line, "$Elements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 1044 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 1045 switch (gmsh->fileFormat) { 1046 case 41: ierr = DMPlexCreateGmsh_ReadElements_v41(gmsh, shift, entities, &numCells, &gmsh_elem);CHKERRQ(ierr); break; 1047 case 40: ierr = DMPlexCreateGmsh_ReadElements_v40(gmsh, shift, entities, &numCells, &gmsh_elem);CHKERRQ(ierr); break; 1048 default: ierr = DMPlexCreateGmsh_ReadElements_v22(gmsh, shift, &numCells, &gmsh_elem);CHKERRQ(ierr); 1049 } 1050 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 1051 ierr = PetscStrncmp(line, "$EndElements", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 1052 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 1053 1054 /* OPTIONAL Read periodic section */ 1055 if (periodic) { 1056 PetscInt pVert, *periodicMapT, *aux; 1057 1058 ierr = PetscMalloc1(numVertices, &periodicMapT);CHKERRQ(ierr); 1059 ierr = PetscBTCreate(numVertices, &periodicV);CHKERRQ(ierr); 1060 for (i = 0; i < numVertices; i++) periodicMapT[i] = i; 1061 1062 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 1063 ierr = PetscStrncmp(line, "$Periodic", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 1064 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 1065 switch (gmsh->fileFormat) { 1066 case 41: ierr = DMPlexCreateGmsh_ReadPeriodic_v41(gmsh, shift, periodicMapT, periodicV);CHKERRQ(ierr); break; 1067 default: ierr = DMPlexCreateGmsh_ReadPeriodic_v40(gmsh, shift, periodicMapT, periodicV);CHKERRQ(ierr); break; 1068 } 1069 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 1070 ierr = PetscStrncmp(line, "$EndPeriodic", PETSC_MAX_PATH_LEN, &match);CHKERRQ(ierr); 1071 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Gmsh file"); 1072 1073 /* we may have slaves of slaves */ 1074 for (i = 0; i < numVertices; i++) { 1075 while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) { 1076 periodicMapT[i] = periodicMapT[periodicMapT[i]]; 1077 } 1078 } 1079 /* periodicMap : from old to new numbering (periodic vertices excluded) 1080 periodicMapI: from new to old numbering */ 1081 ierr = PetscMalloc1(numVertices, &periodicMap);CHKERRQ(ierr); 1082 ierr = PetscMalloc1(numVertices, &periodicMapI);CHKERRQ(ierr); 1083 ierr = PetscMalloc1(numVertices, &aux);CHKERRQ(ierr); 1084 for (i = 0, pVert = 0; i < numVertices; i++) { 1085 if (periodicMapT[i] != i) { 1086 pVert++; 1087 } else { 1088 aux[i] = i - pVert; 1089 periodicMapI[i - pVert] = i; 1090 } 1091 } 1092 for (i = 0 ; i < numVertices; i++) { 1093 periodicMap[i] = aux[periodicMapT[i]]; 1094 } 1095 ierr = PetscFree(periodicMapT);CHKERRQ(ierr); 1096 ierr = PetscFree(aux);CHKERRQ(ierr); 1097 /* remove periodic vertices */ 1098 numVertices = numVertices - pVert; 1099 } 1100 1101 ierr = GmshEntitiesDestroy(&entities);CHKERRQ(ierr); 1102 ierr = PetscFree(gmsh->wbuf);CHKERRQ(ierr); 1103 ierr = PetscFree(gmsh->sbuf);CHKERRQ(ierr); 1104 } 1105 1106 if (parentviewer) { 1107 ierr = PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);CHKERRQ(ierr); 1108 } 1109 1110 if (!rank) { 1111 PetscBool hybrid = PETSC_FALSE; 1112 1113 for (trueNumCells = 0, c = 0; c < numCells; ++c) { 1114 int on = -1; 1115 if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;} 1116 if (gmsh_elem[c].dim == dim) {hybrid = (trueNumCells ? (on != gmsh_elem[c].cellType ? on = gmsh_elem[c].cellType,PETSC_TRUE : hybrid) : (on = gmsh_elem[c].cellType, PETSC_FALSE) ); trueNumCells++;} 1117 /* wedges always indicate an hybrid mesh in PLEX */ 1118 if (on == 6 || on == 13) hybrid = PETSC_TRUE; 1119 } 1120 /* Renumber cells for hybrid grids */ 1121 if (hybrid && enable_hybrid) { 1122 PetscInt hc1 = 0, hc2 = 0, *hybridCells1 = NULL, *hybridCells2 = NULL; 1123 PetscInt cell, tn, *tp; 1124 int n1 = 0,n2 = 0; 1125 1126 ierr = PetscMalloc1(trueNumCells, &hybridCells1);CHKERRQ(ierr); 1127 ierr = PetscMalloc1(trueNumCells, &hybridCells2);CHKERRQ(ierr); 1128 for (cell = 0, c = 0; c < numCells; ++c) { 1129 if (gmsh_elem[c].dim == dim) { 1130 if (!n1) n1 = gmsh_elem[c].cellType; 1131 else if (!n2 && n1 != gmsh_elem[c].cellType) n2 = gmsh_elem[c].cellType; 1132 1133 if (gmsh_elem[c].cellType == n1) hybridCells1[hc1++] = cell; 1134 else if (gmsh_elem[c].cellType == n2) hybridCells2[hc2++] = cell; 1135 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle more than 2 cell types"); 1136 cell++; 1137 } 1138 } 1139 1140 switch (n1) { 1141 case 2: /* triangles */ 1142 case 9: 1143 switch (n2) { 1144 case 0: /* single type mesh */ 1145 case 3: /* quads */ 1146 case 10: 1147 break; 1148 default: 1149 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1150 } 1151 break; 1152 case 3: 1153 case 10: 1154 switch (n2) { 1155 case 0: /* single type mesh */ 1156 case 2: /* swap since we list simplices first */ 1157 case 9: 1158 tn = hc1; 1159 hc1 = hc2; 1160 hc2 = tn; 1161 1162 tp = hybridCells1; 1163 hybridCells1 = hybridCells2; 1164 hybridCells2 = tp; 1165 break; 1166 default: 1167 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1168 } 1169 break; 1170 case 4: /* tetrahedra */ 1171 case 11: 1172 switch (n2) { 1173 case 0: /* single type mesh */ 1174 case 6: /* wedges */ 1175 case 13: 1176 break; 1177 default: 1178 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1179 } 1180 break; 1181 case 6: 1182 case 13: 1183 switch (n2) { 1184 case 0: /* single type mesh */ 1185 case 4: /* swap since we list simplices first */ 1186 case 11: 1187 tn = hc1; 1188 hc1 = hc2; 1189 hc2 = tn; 1190 1191 tp = hybridCells1; 1192 hybridCells1 = hybridCells2; 1193 hybridCells2 = tp; 1194 break; 1195 default: 1196 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1197 } 1198 break; 1199 default: 1200 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1201 } 1202 cMax = hc1; 1203 ierr = PetscMalloc1(trueNumCells, &hybridMap);CHKERRQ(ierr); 1204 for (cell = 0; cell < hc1; cell++) hybridMap[hybridCells1[cell]] = cell; 1205 for (cell = 0; cell < hc2; cell++) hybridMap[hybridCells2[cell]] = cell + hc1; 1206 ierr = PetscFree(hybridCells1);CHKERRQ(ierr); 1207 ierr = PetscFree(hybridCells2);CHKERRQ(ierr); 1208 } 1209 1210 } 1211 1212 /* Allocate the cell-vertex mesh */ 1213 ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr); 1214 for (cell = 0, c = 0; c < numCells; ++c) { 1215 if (gmsh_elem[c].dim == dim) { 1216 ierr = DMPlexSetConeSize(*dm, hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].numNodes);CHKERRQ(ierr); 1217 cell++; 1218 } 1219 } 1220 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1221 /* Add cell-vertex connections */ 1222 for (cell = 0, c = 0; c < numCells; ++c) { 1223 if (gmsh_elem[c].dim == dim) { 1224 PetscInt pcone[8], corner; 1225 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1226 const PetscInt cc = gmsh_elem[c].nodes[corner] - shift; 1227 pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells; 1228 } 1229 if (dim == 3) { 1230 /* Tetrahedra are inverted */ 1231 if (gmsh_elem[c].cellType == 4) { 1232 PetscInt tmp = pcone[0]; 1233 pcone[0] = pcone[1]; 1234 pcone[1] = tmp; 1235 } 1236 /* Hexahedra are inverted */ 1237 if (gmsh_elem[c].cellType == 5) { 1238 PetscInt tmp = pcone[1]; 1239 pcone[1] = pcone[3]; 1240 pcone[3] = tmp; 1241 } 1242 /* Prisms are inverted */ 1243 if (gmsh_elem[c].cellType == 6) { 1244 PetscInt tmp; 1245 1246 tmp = pcone[1]; 1247 pcone[1] = pcone[2]; 1248 pcone[2] = tmp; 1249 tmp = pcone[4]; 1250 pcone[4] = pcone[5]; 1251 pcone[5] = tmp; 1252 } 1253 } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */ 1254 /* quads are input to PLEX as prisms */ 1255 if (gmsh_elem[c].cellType == 3) { 1256 PetscInt tmp = pcone[2]; 1257 pcone[2] = pcone[3]; 1258 pcone[3] = tmp; 1259 } 1260 } 1261 ierr = DMPlexSetCone(*dm, hybridMap ? hybridMap[cell] : cell, pcone);CHKERRQ(ierr); 1262 cell++; 1263 } 1264 } 1265 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 1266 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1267 ierr = DMPlexSetHybridBounds(*dm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1268 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1269 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1270 if (interpolate) { 1271 DM idm; 1272 1273 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1274 ierr = DMDestroy(dm);CHKERRQ(ierr); 1275 *dm = idm; 1276 } 1277 1278 if (usemarker && !interpolate && dim > 1) SETERRQ(comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation"); 1279 if (!rank && usemarker) { 1280 PetscInt f, fStart, fEnd; 1281 1282 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 1283 ierr = DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1284 for (f = fStart; f < fEnd; ++f) { 1285 PetscInt suppSize; 1286 1287 ierr = DMPlexGetSupportSize(*dm, f, &suppSize);CHKERRQ(ierr); 1288 if (suppSize == 1) { 1289 PetscInt *cone = NULL, coneSize, p; 1290 1291 ierr = DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);CHKERRQ(ierr); 1292 for (p = 0; p < coneSize; p += 2) { 1293 ierr = DMSetLabelValue(*dm, "marker", cone[p], 1);CHKERRQ(ierr); 1294 } 1295 ierr = DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);CHKERRQ(ierr); 1296 } 1297 } 1298 } 1299 1300 if (!rank) { 1301 PetscInt vStart, vEnd; 1302 1303 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1304 for (cell = 0, c = 0; c < numCells; ++c) { 1305 1306 /* Create face sets */ 1307 if (interpolate && gmsh_elem[c].dim == dim-1) { 1308 const PetscInt *join; 1309 PetscInt joinSize, pcone[8], corner; 1310 /* Find the relevant facet with vertex joins */ 1311 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1312 const PetscInt cc = gmsh_elem[c].nodes[corner] - shift; 1313 pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart; 1314 } 1315 ierr = DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, pcone, &joinSize, &join);CHKERRQ(ierr); 1316 if (joinSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for GMsh element %d (Plex cell %D)", gmsh_elem[c].id, c); 1317 ierr = DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);CHKERRQ(ierr); 1318 ierr = DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 1319 } 1320 1321 /* Create cell sets */ 1322 if (gmsh_elem[c].dim == dim) { 1323 if (gmsh_elem[c].numTags > 0) { 1324 ierr = DMSetLabelValue(*dm, "Cell Sets", hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].tags[0]);CHKERRQ(ierr); 1325 } 1326 cell++; 1327 } 1328 1329 /* Create vertex sets */ 1330 if (gmsh_elem[c].dim == 0) { 1331 if (gmsh_elem[c].numTags > 0) { 1332 const PetscInt cc = gmsh_elem[c].nodes[0] - shift; 1333 const PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart; 1334 ierr = DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);CHKERRQ(ierr); 1335 } 1336 } 1337 } 1338 } 1339 1340 /* Create coordinates */ 1341 if (embedDim < 0) embedDim = dim; 1342 ierr = DMSetCoordinateDim(*dm, embedDim);CHKERRQ(ierr); 1343 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1344 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1345 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 1346 if (periodicMap) { /* we need to localize coordinates on cells */ 1347 ierr = PetscSectionSetChart(coordSection, 0, trueNumCells + numVertices);CHKERRQ(ierr); 1348 } else { 1349 ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr); 1350 } 1351 for (v = trueNumCells; v < trueNumCells+numVertices; ++v) { 1352 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 1353 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 1354 } 1355 if (periodicMap) { 1356 ierr = PetscBTCreate(trueNumCells, &periodicC);CHKERRQ(ierr); 1357 for (cell = 0, c = 0; c < numCells; ++c) { 1358 if (gmsh_elem[c].dim == dim) { 1359 PetscInt corner; 1360 PetscBool pc = PETSC_FALSE; 1361 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1362 pc = (PetscBool)(pc || PetscBTLookup(periodicV, gmsh_elem[c].nodes[corner] - shift)); 1363 } 1364 if (pc) { 1365 PetscInt dof = gmsh_elem[c].numNodes*embedDim; 1366 PetscInt ucell = hybridMap ? hybridMap[cell] : cell; 1367 ierr = PetscBTSet(periodicC, ucell);CHKERRQ(ierr); 1368 ierr = PetscSectionSetDof(coordSection, ucell, dof);CHKERRQ(ierr); 1369 ierr = PetscSectionSetFieldDof(coordSection, ucell, 0, dof);CHKERRQ(ierr); 1370 } 1371 cell++; 1372 } 1373 } 1374 } 1375 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1376 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1377 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1378 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1379 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1380 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 1381 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 1382 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1383 if (periodicMap) { 1384 PetscInt off; 1385 1386 for (cell = 0, c = 0; c < numCells; ++c) { 1387 PetscInt pcone[8], corner; 1388 if (gmsh_elem[c].dim == dim) { 1389 PetscInt ucell = hybridMap ? hybridMap[cell] : cell; 1390 if (PetscUnlikely(PetscBTLookup(periodicC, ucell))) { 1391 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1392 pcone[corner] = gmsh_elem[c].nodes[corner] - shift; 1393 } 1394 if (dim == 3) { 1395 /* Tetrahedra are inverted */ 1396 if (gmsh_elem[c].cellType == 4) { 1397 PetscInt tmp = pcone[0]; 1398 pcone[0] = pcone[1]; 1399 pcone[1] = tmp; 1400 } 1401 /* Hexahedra are inverted */ 1402 if (gmsh_elem[c].cellType == 5) { 1403 PetscInt tmp = pcone[1]; 1404 pcone[1] = pcone[3]; 1405 pcone[3] = tmp; 1406 } 1407 /* Prisms are inverted */ 1408 if (gmsh_elem[c].cellType == 6) { 1409 PetscInt tmp; 1410 1411 tmp = pcone[1]; 1412 pcone[1] = pcone[2]; 1413 pcone[2] = tmp; 1414 tmp = pcone[4]; 1415 pcone[4] = pcone[5]; 1416 pcone[5] = tmp; 1417 } 1418 } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */ 1419 /* quads are input to PLEX as prisms */ 1420 if (gmsh_elem[c].cellType == 3) { 1421 PetscInt tmp = pcone[2]; 1422 pcone[2] = pcone[3]; 1423 pcone[3] = tmp; 1424 } 1425 } 1426 ierr = PetscSectionGetOffset(coordSection, ucell, &off);CHKERRQ(ierr); 1427 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1428 v = pcone[corner]; 1429 for (d = 0; d < embedDim; ++d) { 1430 coords[off++] = (PetscReal) coordsIn[v*3+d]; 1431 } 1432 } 1433 } 1434 cell++; 1435 } 1436 } 1437 for (v = 0; v < numVertices; ++v) { 1438 ierr = PetscSectionGetOffset(coordSection, v + trueNumCells, &off);CHKERRQ(ierr); 1439 for (d = 0; d < embedDim; ++d) { 1440 coords[off+d] = (PetscReal) coordsIn[periodicMapI[v]*3+d]; 1441 } 1442 } 1443 } else { 1444 for (v = 0; v < numVertices; ++v) { 1445 for (d = 0; d < embedDim; ++d) { 1446 coords[v*embedDim+d] = (PetscReal) coordsIn[v*3+d]; 1447 } 1448 } 1449 } 1450 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1451 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1452 ierr = DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);CHKERRQ(ierr); 1453 1454 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 1455 ierr = PetscFree(gmsh_elem);CHKERRQ(ierr); 1456 ierr = PetscFree(hybridMap);CHKERRQ(ierr); 1457 ierr = PetscFree(periodicMap);CHKERRQ(ierr); 1458 ierr = PetscFree(periodicMapI);CHKERRQ(ierr); 1459 ierr = PetscBTDestroy(&periodicV);CHKERRQ(ierr); 1460 ierr = PetscBTDestroy(&periodicC);CHKERRQ(ierr); 1461 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1462 1463 ierr = PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 1464 PetscFunctionReturn(0); 1465 } 1466