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 .seealso: DMPLEX, DMCreate() 1010 @*/ 1011 PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm) 1012 { 1013 PetscViewer parentviewer = NULL; 1014 double *coordsIn = NULL; 1015 GmshEntities *entities = NULL; 1016 GmshElement *gmsh_elem = NULL; 1017 PetscSection coordSection; 1018 Vec coordinates; 1019 PetscBT periodicV = NULL, periodicC = NULL; 1020 PetscScalar *coords; 1021 PetscInt dim = 0, embedDim = -1, coordSize, c, v, d, cell, *periodicMap = NULL, *periodicMapI = NULL, *hybridMap = NULL, cMax = PETSC_DETERMINE; 1022 PetscInt numVertices = 0, numCells = 0, trueNumCells = 0; 1023 int i, shift = 1; 1024 PetscMPIInt rank; 1025 PetscBool binary, zerobase = PETSC_FALSE, usemarker = PETSC_FALSE; 1026 PetscBool enable_hybrid = interpolate, periodic = PETSC_TRUE; 1027 PetscErrorCode ierr; 1028 1029 PetscFunctionBegin; 1030 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1031 ierr = PetscObjectOptionsBegin((PetscObject)viewer);CHKERRQ(ierr); 1032 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Gmsh options");CHKERRQ(ierr); 1033 ierr = PetscOptionsBool("-dm_plex_gmsh_hybrid", "Generate hybrid cell bounds", "DMPlexCreateGmsh", enable_hybrid, &enable_hybrid, NULL);CHKERRQ(ierr); 1034 ierr = PetscOptionsBool("-dm_plex_gmsh_periodic","Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL);CHKERRQ(ierr); 1035 ierr = PetscOptionsBool("-dm_plex_gmsh_use_marker", "Generate marker label", "DMPlexCreateGmsh", usemarker, &usemarker, NULL);CHKERRQ(ierr); 1036 ierr = PetscOptionsBool("-dm_plex_gmsh_zero_base", "Read Gmsh file with zero base indices", "DMPlexCreateGmsh", zerobase, &zerobase, NULL);CHKERRQ(ierr); 1037 ierr = PetscOptionsInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", embedDim, &embedDim, NULL);CHKERRQ(ierr); 1038 ierr = PetscOptionsTail();CHKERRQ(ierr); 1039 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1040 if (zerobase) shift = 0; 1041 1042 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1043 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1044 ierr = PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 1045 1046 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);CHKERRQ(ierr); 1047 1048 /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */ 1049 if (binary) { 1050 parentviewer = viewer; 1051 ierr = PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);CHKERRQ(ierr); 1052 } 1053 1054 if (!rank) { 1055 GmshFile gmsh_, *gmsh = &gmsh_; 1056 char line[PETSC_MAX_PATH_LEN]; 1057 PetscBool match; 1058 1059 ierr = PetscMemzero(gmsh,sizeof(GmshFile));CHKERRQ(ierr); 1060 gmsh->viewer = viewer; 1061 gmsh->binary = binary; 1062 1063 /* Read mesh format */ 1064 ierr = GmshReadSection(gmsh, line);CHKERRQ(ierr); 1065 ierr = GmshExpect(gmsh, "$MeshFormat", line);CHKERRQ(ierr); 1066 ierr = DMPlexCreateGmsh_ReadMeshFormat(gmsh);CHKERRQ(ierr); 1067 ierr = GmshReadEndSection(gmsh, "$EndMeshFormat", line);CHKERRQ(ierr); 1068 1069 /* OPTIONAL Read physical names */ 1070 ierr = GmshReadSection(gmsh, line);CHKERRQ(ierr); 1071 ierr = GmshMatch(gmsh,"$PhysicalNames", line, &match);CHKERRQ(ierr); 1072 if (match) { 1073 ierr = DMPlexCreateGmsh_ReadPhysicalNames(gmsh);CHKERRQ(ierr); 1074 ierr = GmshReadEndSection(gmsh, "$EndPhysicalNames", line);CHKERRQ(ierr); 1075 /* Initial read for entity section */ 1076 ierr = GmshReadSection(gmsh, line);CHKERRQ(ierr); 1077 } 1078 1079 /* Read entities */ 1080 if (gmsh->fileFormat >= 40) { 1081 ierr = GmshExpect(gmsh, "$Entities", line);CHKERRQ(ierr); 1082 switch (gmsh->fileFormat) { 1083 case 41: ierr = DMPlexCreateGmsh_ReadEntities_v41(gmsh, &entities);CHKERRQ(ierr); break; 1084 default: ierr = DMPlexCreateGmsh_ReadEntities_v40(gmsh, &entities);CHKERRQ(ierr); break; 1085 } 1086 ierr = GmshReadEndSection(gmsh, "$EndEntities", line);CHKERRQ(ierr); 1087 /* Initial read for nodes section */ 1088 ierr = GmshReadSection(gmsh, line);CHKERRQ(ierr); 1089 } 1090 1091 /* Read nodes */ 1092 ierr = GmshExpect(gmsh, "$Nodes", line);CHKERRQ(ierr); 1093 switch (gmsh->fileFormat) { 1094 case 41: ierr = DMPlexCreateGmsh_ReadNodes_v41(gmsh, shift, &numVertices, &coordsIn);CHKERRQ(ierr); break; 1095 case 40: ierr = DMPlexCreateGmsh_ReadNodes_v40(gmsh, shift, &numVertices, &coordsIn);CHKERRQ(ierr); break; 1096 default: ierr = DMPlexCreateGmsh_ReadNodes_v22(gmsh, shift, &numVertices, &coordsIn);CHKERRQ(ierr); break; 1097 } 1098 ierr = GmshReadEndSection(gmsh, "$EndNodes", line);CHKERRQ(ierr); 1099 1100 /* Read elements */ 1101 ierr = GmshReadSection(gmsh, line);CHKERRQ(ierr);; 1102 ierr = GmshExpect(gmsh, "$Elements", line);CHKERRQ(ierr); 1103 switch (gmsh->fileFormat) { 1104 case 41: ierr = DMPlexCreateGmsh_ReadElements_v41(gmsh, shift, entities, &numCells, &gmsh_elem);CHKERRQ(ierr); break; 1105 case 40: ierr = DMPlexCreateGmsh_ReadElements_v40(gmsh, shift, entities, &numCells, &gmsh_elem);CHKERRQ(ierr); break; 1106 default: ierr = DMPlexCreateGmsh_ReadElements_v22(gmsh, shift, &numCells, &gmsh_elem);CHKERRQ(ierr); break; 1107 } 1108 ierr = GmshReadEndSection(gmsh, "$EndElements", line);CHKERRQ(ierr); 1109 1110 /* OPTIONAL Read periodic section */ 1111 if (periodic) { 1112 ierr = GmshReadSection(gmsh, line);CHKERRQ(ierr); 1113 ierr = GmshMatch(gmsh, "$Periodic", line, &periodic);CHKERRQ(ierr); 1114 } 1115 if (periodic) { 1116 PetscInt pVert, *periodicMapT, *aux; 1117 1118 ierr = PetscMalloc1(numVertices, &periodicMapT);CHKERRQ(ierr); 1119 ierr = PetscBTCreate(numVertices, &periodicV);CHKERRQ(ierr); 1120 for (i = 0; i < numVertices; i++) periodicMapT[i] = i; 1121 1122 ierr = GmshExpect(gmsh, "$Periodic", line);CHKERRQ(ierr); 1123 switch (gmsh->fileFormat) { 1124 case 41: ierr = DMPlexCreateGmsh_ReadPeriodic_v41(gmsh, shift, periodicMapT, periodicV);CHKERRQ(ierr); break; 1125 default: ierr = DMPlexCreateGmsh_ReadPeriodic_v40(gmsh, shift, periodicMapT, periodicV);CHKERRQ(ierr); break; 1126 } 1127 ierr = GmshReadEndSection(gmsh, "$EndPeriodic", line);CHKERRQ(ierr); 1128 1129 /* we may have slaves of slaves */ 1130 for (i = 0; i < numVertices; i++) { 1131 while (periodicMapT[periodicMapT[i]] != periodicMapT[i]) { 1132 periodicMapT[i] = periodicMapT[periodicMapT[i]]; 1133 } 1134 } 1135 /* periodicMap : from old to new numbering (periodic vertices excluded) 1136 periodicMapI: from new to old numbering */ 1137 ierr = PetscMalloc1(numVertices, &periodicMap);CHKERRQ(ierr); 1138 ierr = PetscMalloc1(numVertices, &periodicMapI);CHKERRQ(ierr); 1139 ierr = PetscMalloc1(numVertices, &aux);CHKERRQ(ierr); 1140 for (i = 0, pVert = 0; i < numVertices; i++) { 1141 if (periodicMapT[i] != i) { 1142 pVert++; 1143 } else { 1144 aux[i] = i - pVert; 1145 periodicMapI[i - pVert] = i; 1146 } 1147 } 1148 for (i = 0 ; i < numVertices; i++) { 1149 periodicMap[i] = aux[periodicMapT[i]]; 1150 } 1151 ierr = PetscFree(periodicMapT);CHKERRQ(ierr); 1152 ierr = PetscFree(aux);CHKERRQ(ierr); 1153 /* remove periodic vertices */ 1154 numVertices = numVertices - pVert; 1155 } 1156 1157 ierr = GmshEntitiesDestroy(&entities);CHKERRQ(ierr); 1158 ierr = PetscFree(gmsh->wbuf);CHKERRQ(ierr); 1159 ierr = PetscFree(gmsh->sbuf);CHKERRQ(ierr); 1160 } 1161 1162 if (parentviewer) { 1163 ierr = PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);CHKERRQ(ierr); 1164 } 1165 1166 if (!rank) { 1167 PetscBool hybrid = PETSC_FALSE; 1168 PetscInt cellType = -1; 1169 1170 for (trueNumCells = 0, c = 0; c < numCells; ++c) { 1171 if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0; hybrid = PETSC_FALSE; cellType = -1;} 1172 if (gmsh_elem[c].dim < dim) continue; 1173 if (cellType == -1) cellType = gmsh_elem[c].cellType; 1174 /* different cell type indicate an hybrid mesh in PLEX */ 1175 if (cellType != gmsh_elem[c].cellType) hybrid = PETSC_TRUE; 1176 /* wedges always indicate an hybrid mesh in PLEX */ 1177 if (cellType == 6 || cellType == 13) hybrid = PETSC_TRUE; 1178 trueNumCells++; 1179 } 1180 /* Renumber cells for hybrid grids */ 1181 if (hybrid && enable_hybrid) { 1182 PetscInt hc1 = 0, hc2 = 0, *hybridCells1 = NULL, *hybridCells2 = NULL; 1183 PetscInt cell, tn, *tp; 1184 int n1 = 0,n2 = 0; 1185 1186 ierr = PetscMalloc1(trueNumCells, &hybridCells1);CHKERRQ(ierr); 1187 ierr = PetscMalloc1(trueNumCells, &hybridCells2);CHKERRQ(ierr); 1188 for (cell = 0, c = 0; c < numCells; ++c) { 1189 if (gmsh_elem[c].dim == dim) { 1190 if (!n1) n1 = gmsh_elem[c].cellType; 1191 else if (!n2 && n1 != gmsh_elem[c].cellType) n2 = gmsh_elem[c].cellType; 1192 1193 if (gmsh_elem[c].cellType == n1) hybridCells1[hc1++] = cell; 1194 else if (gmsh_elem[c].cellType == n2) hybridCells2[hc2++] = cell; 1195 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle more than 2 cell types"); 1196 cell++; 1197 } 1198 } 1199 1200 switch (n1) { 1201 case 2: /* triangles */ 1202 case 9: 1203 switch (n2) { 1204 case 0: /* single type mesh */ 1205 case 3: /* quads */ 1206 case 10: 1207 break; 1208 default: 1209 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1210 } 1211 break; 1212 case 3: 1213 case 10: 1214 switch (n2) { 1215 case 0: /* single type mesh */ 1216 case 2: /* swap since we list simplices first */ 1217 case 9: 1218 tn = hc1; 1219 hc1 = hc2; 1220 hc2 = tn; 1221 1222 tp = hybridCells1; 1223 hybridCells1 = hybridCells2; 1224 hybridCells2 = tp; 1225 break; 1226 default: 1227 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1228 } 1229 break; 1230 case 4: /* tetrahedra */ 1231 case 11: 1232 switch (n2) { 1233 case 0: /* single type mesh */ 1234 case 6: /* wedges */ 1235 case 13: 1236 break; 1237 default: 1238 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1239 } 1240 break; 1241 case 6: 1242 case 13: 1243 switch (n2) { 1244 case 0: /* single type mesh */ 1245 case 4: /* swap since we list simplices first */ 1246 case 11: 1247 tn = hc1; 1248 hc1 = hc2; 1249 hc2 = tn; 1250 1251 tp = hybridCells1; 1252 hybridCells1 = hybridCells2; 1253 hybridCells2 = tp; 1254 break; 1255 default: 1256 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1257 } 1258 break; 1259 default: 1260 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1261 } 1262 cMax = hc1; 1263 ierr = PetscMalloc1(trueNumCells, &hybridMap);CHKERRQ(ierr); 1264 for (cell = 0; cell < hc1; cell++) hybridMap[hybridCells1[cell]] = cell; 1265 for (cell = 0; cell < hc2; cell++) hybridMap[hybridCells2[cell]] = cell + hc1; 1266 ierr = PetscFree(hybridCells1);CHKERRQ(ierr); 1267 ierr = PetscFree(hybridCells2);CHKERRQ(ierr); 1268 } 1269 1270 } 1271 1272 /* Allocate the cell-vertex mesh */ 1273 ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr); 1274 for (cell = 0, c = 0; c < numCells; ++c) { 1275 if (gmsh_elem[c].dim == dim) { 1276 ierr = DMPlexSetConeSize(*dm, hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].numNodes);CHKERRQ(ierr); 1277 cell++; 1278 } 1279 } 1280 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1281 /* Add cell-vertex connections */ 1282 for (cell = 0, c = 0; c < numCells; ++c) { 1283 if (gmsh_elem[c].dim == dim) { 1284 PetscInt pcone[8], corner; 1285 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1286 const PetscInt cc = gmsh_elem[c].nodes[corner] - shift; 1287 pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells; 1288 } 1289 if (dim == 3) { 1290 /* Tetrahedra are inverted */ 1291 if (gmsh_elem[c].cellType == 4) { 1292 PetscInt tmp = pcone[0]; 1293 pcone[0] = pcone[1]; 1294 pcone[1] = tmp; 1295 } 1296 /* Hexahedra are inverted */ 1297 if (gmsh_elem[c].cellType == 5) { 1298 PetscInt tmp = pcone[1]; 1299 pcone[1] = pcone[3]; 1300 pcone[3] = tmp; 1301 } 1302 /* Prisms are inverted */ 1303 if (gmsh_elem[c].cellType == 6) { 1304 PetscInt tmp; 1305 1306 tmp = pcone[1]; 1307 pcone[1] = pcone[2]; 1308 pcone[2] = tmp; 1309 tmp = pcone[4]; 1310 pcone[4] = pcone[5]; 1311 pcone[5] = tmp; 1312 } 1313 } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */ 1314 /* quads are input to PLEX as prisms */ 1315 if (gmsh_elem[c].cellType == 3) { 1316 PetscInt tmp = pcone[2]; 1317 pcone[2] = pcone[3]; 1318 pcone[3] = tmp; 1319 } 1320 } 1321 ierr = DMPlexSetCone(*dm, hybridMap ? hybridMap[cell] : cell, pcone);CHKERRQ(ierr); 1322 cell++; 1323 } 1324 } 1325 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 1326 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1327 ierr = DMPlexSetHybridBounds(*dm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1328 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1329 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1330 if (interpolate) { 1331 DM idm; 1332 1333 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1334 ierr = DMDestroy(dm);CHKERRQ(ierr); 1335 *dm = idm; 1336 } 1337 1338 if (usemarker && !interpolate && dim > 1) SETERRQ(comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation"); 1339 if (!rank && usemarker) { 1340 PetscInt f, fStart, fEnd; 1341 1342 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 1343 ierr = DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1344 for (f = fStart; f < fEnd; ++f) { 1345 PetscInt suppSize; 1346 1347 ierr = DMPlexGetSupportSize(*dm, f, &suppSize);CHKERRQ(ierr); 1348 if (suppSize == 1) { 1349 PetscInt *cone = NULL, coneSize, p; 1350 1351 ierr = DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);CHKERRQ(ierr); 1352 for (p = 0; p < coneSize; p += 2) { 1353 ierr = DMSetLabelValue(*dm, "marker", cone[p], 1);CHKERRQ(ierr); 1354 } 1355 ierr = DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);CHKERRQ(ierr); 1356 } 1357 } 1358 } 1359 1360 if (!rank) { 1361 PetscInt vStart, vEnd; 1362 1363 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1364 for (cell = 0, c = 0; c < numCells; ++c) { 1365 1366 /* Create face sets */ 1367 if (interpolate && gmsh_elem[c].dim == dim-1) { 1368 const PetscInt *join; 1369 PetscInt joinSize, pcone[8], corner; 1370 /* Find the relevant facet with vertex joins */ 1371 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1372 const PetscInt cc = gmsh_elem[c].nodes[corner] - shift; 1373 pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart; 1374 } 1375 ierr = DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, pcone, &joinSize, &join);CHKERRQ(ierr); 1376 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); 1377 ierr = DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);CHKERRQ(ierr); 1378 ierr = DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 1379 } 1380 1381 /* Create cell sets */ 1382 if (gmsh_elem[c].dim == dim) { 1383 if (gmsh_elem[c].numTags > 0) { 1384 ierr = DMSetLabelValue(*dm, "Cell Sets", hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].tags[0]);CHKERRQ(ierr); 1385 } 1386 cell++; 1387 } 1388 1389 /* Create vertex sets */ 1390 if (gmsh_elem[c].dim == 0) { 1391 if (gmsh_elem[c].numTags > 0) { 1392 const PetscInt cc = gmsh_elem[c].nodes[0] - shift; 1393 const PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart; 1394 ierr = DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);CHKERRQ(ierr); 1395 } 1396 } 1397 } 1398 } 1399 1400 /* Create coordinates */ 1401 if (embedDim < 0) embedDim = dim; 1402 ierr = DMSetCoordinateDim(*dm, embedDim);CHKERRQ(ierr); 1403 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1404 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1405 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 1406 if (periodicMap) { /* we need to localize coordinates on cells */ 1407 ierr = PetscSectionSetChart(coordSection, 0, trueNumCells + numVertices);CHKERRQ(ierr); 1408 } else { 1409 ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr); 1410 } 1411 for (v = trueNumCells; v < trueNumCells+numVertices; ++v) { 1412 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 1413 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 1414 } 1415 if (periodicMap) { 1416 ierr = PetscBTCreate(trueNumCells, &periodicC);CHKERRQ(ierr); 1417 for (cell = 0, c = 0; c < numCells; ++c) { 1418 if (gmsh_elem[c].dim == dim) { 1419 PetscInt corner; 1420 PetscBool pc = PETSC_FALSE; 1421 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1422 pc = (PetscBool)(pc || PetscBTLookup(periodicV, gmsh_elem[c].nodes[corner] - shift)); 1423 } 1424 if (pc) { 1425 PetscInt dof = gmsh_elem[c].numNodes*embedDim; 1426 PetscInt ucell = hybridMap ? hybridMap[cell] : cell; 1427 ierr = PetscBTSet(periodicC, ucell);CHKERRQ(ierr); 1428 ierr = PetscSectionSetDof(coordSection, ucell, dof);CHKERRQ(ierr); 1429 ierr = PetscSectionSetFieldDof(coordSection, ucell, 0, dof);CHKERRQ(ierr); 1430 } 1431 cell++; 1432 } 1433 } 1434 } 1435 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1436 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1437 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1438 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1439 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1440 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 1441 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 1442 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1443 if (periodicMap) { 1444 PetscInt off; 1445 1446 for (cell = 0, c = 0; c < numCells; ++c) { 1447 PetscInt pcone[8], corner; 1448 if (gmsh_elem[c].dim == dim) { 1449 PetscInt ucell = hybridMap ? hybridMap[cell] : cell; 1450 if (PetscUnlikely(PetscBTLookup(periodicC, ucell))) { 1451 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1452 pcone[corner] = gmsh_elem[c].nodes[corner] - shift; 1453 } 1454 if (dim == 3) { 1455 /* Tetrahedra are inverted */ 1456 if (gmsh_elem[c].cellType == 4) { 1457 PetscInt tmp = pcone[0]; 1458 pcone[0] = pcone[1]; 1459 pcone[1] = tmp; 1460 } 1461 /* Hexahedra are inverted */ 1462 if (gmsh_elem[c].cellType == 5) { 1463 PetscInt tmp = pcone[1]; 1464 pcone[1] = pcone[3]; 1465 pcone[3] = tmp; 1466 } 1467 /* Prisms are inverted */ 1468 if (gmsh_elem[c].cellType == 6) { 1469 PetscInt tmp; 1470 1471 tmp = pcone[1]; 1472 pcone[1] = pcone[2]; 1473 pcone[2] = tmp; 1474 tmp = pcone[4]; 1475 pcone[4] = pcone[5]; 1476 pcone[5] = tmp; 1477 } 1478 } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */ 1479 /* quads are input to PLEX as prisms */ 1480 if (gmsh_elem[c].cellType == 3) { 1481 PetscInt tmp = pcone[2]; 1482 pcone[2] = pcone[3]; 1483 pcone[3] = tmp; 1484 } 1485 } 1486 ierr = PetscSectionGetOffset(coordSection, ucell, &off);CHKERRQ(ierr); 1487 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1488 v = pcone[corner]; 1489 for (d = 0; d < embedDim; ++d) { 1490 coords[off++] = (PetscReal) coordsIn[v*3+d]; 1491 } 1492 } 1493 } 1494 cell++; 1495 } 1496 } 1497 for (v = 0; v < numVertices; ++v) { 1498 ierr = PetscSectionGetOffset(coordSection, v + trueNumCells, &off);CHKERRQ(ierr); 1499 for (d = 0; d < embedDim; ++d) { 1500 coords[off+d] = (PetscReal) coordsIn[periodicMapI[v]*3+d]; 1501 } 1502 } 1503 } else { 1504 for (v = 0; v < numVertices; ++v) { 1505 for (d = 0; d < embedDim; ++d) { 1506 coords[v*embedDim+d] = (PetscReal) coordsIn[v*3+d]; 1507 } 1508 } 1509 } 1510 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1511 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1512 1513 periodic = periodicMap ? PETSC_TRUE : PETSC_FALSE; 1514 ierr = MPI_Bcast(&periodic, 1, MPIU_BOOL, 0, comm);CHKERRQ(ierr); 1515 ierr = DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);CHKERRQ(ierr); 1516 1517 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 1518 ierr = PetscFree(gmsh_elem);CHKERRQ(ierr); 1519 ierr = PetscFree(hybridMap);CHKERRQ(ierr); 1520 ierr = PetscFree(periodicMap);CHKERRQ(ierr); 1521 ierr = PetscFree(periodicMapI);CHKERRQ(ierr); 1522 ierr = PetscBTDestroy(&periodicV);CHKERRQ(ierr); 1523 ierr = PetscBTDestroy(&periodicC);CHKERRQ(ierr); 1524 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1525 1526 ierr = PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 1527 PetscFunctionReturn(0); 1528 } 1529