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