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