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