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