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