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 = PETSC_FALSE, 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 1172 for (trueNumCells = 0, c = 0; c < numCells; ++c) { 1173 int on = -1; 1174 if (gmsh_elem[c].dim > dim) {dim = gmsh_elem[c].dim; trueNumCells = 0;} 1175 if (gmsh_elem[c].dim == dim) {hybrid = (trueNumCells ? (on != gmsh_elem[c].cellType ? on = gmsh_elem[c].cellType,PETSC_TRUE : hybrid) : (on = gmsh_elem[c].cellType, PETSC_FALSE) ); trueNumCells++;} 1176 /* wedges always indicate an hybrid mesh in PLEX */ 1177 if (on == 6 || on == 13) hybrid = PETSC_TRUE; 1178 } 1179 /* Renumber cells for hybrid grids */ 1180 if (hybrid && enable_hybrid) { 1181 PetscInt hc1 = 0, hc2 = 0, *hybridCells1 = NULL, *hybridCells2 = NULL; 1182 PetscInt cell, tn, *tp; 1183 int n1 = 0,n2 = 0; 1184 1185 ierr = PetscMalloc1(trueNumCells, &hybridCells1);CHKERRQ(ierr); 1186 ierr = PetscMalloc1(trueNumCells, &hybridCells2);CHKERRQ(ierr); 1187 for (cell = 0, c = 0; c < numCells; ++c) { 1188 if (gmsh_elem[c].dim == dim) { 1189 if (!n1) n1 = gmsh_elem[c].cellType; 1190 else if (!n2 && n1 != gmsh_elem[c].cellType) n2 = gmsh_elem[c].cellType; 1191 1192 if (gmsh_elem[c].cellType == n1) hybridCells1[hc1++] = cell; 1193 else if (gmsh_elem[c].cellType == n2) hybridCells2[hc2++] = cell; 1194 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle more than 2 cell types"); 1195 cell++; 1196 } 1197 } 1198 1199 switch (n1) { 1200 case 2: /* triangles */ 1201 case 9: 1202 switch (n2) { 1203 case 0: /* single type mesh */ 1204 case 3: /* quads */ 1205 case 10: 1206 break; 1207 default: 1208 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1209 } 1210 break; 1211 case 3: 1212 case 10: 1213 switch (n2) { 1214 case 0: /* single type mesh */ 1215 case 2: /* swap since we list simplices first */ 1216 case 9: 1217 tn = hc1; 1218 hc1 = hc2; 1219 hc2 = tn; 1220 1221 tp = hybridCells1; 1222 hybridCells1 = hybridCells2; 1223 hybridCells2 = tp; 1224 break; 1225 default: 1226 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1227 } 1228 break; 1229 case 4: /* tetrahedra */ 1230 case 11: 1231 switch (n2) { 1232 case 0: /* single type mesh */ 1233 case 6: /* wedges */ 1234 case 13: 1235 break; 1236 default: 1237 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1238 } 1239 break; 1240 case 6: 1241 case 13: 1242 switch (n2) { 1243 case 0: /* single type mesh */ 1244 case 4: /* swap since we list simplices first */ 1245 case 11: 1246 tn = hc1; 1247 hc1 = hc2; 1248 hc2 = tn; 1249 1250 tp = hybridCells1; 1251 hybridCells1 = hybridCells2; 1252 hybridCells2 = tp; 1253 break; 1254 default: 1255 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1256 } 1257 break; 1258 default: 1259 SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported cell types %d and %d",n1, n2); 1260 } 1261 cMax = hc1; 1262 ierr = PetscMalloc1(trueNumCells, &hybridMap);CHKERRQ(ierr); 1263 for (cell = 0; cell < hc1; cell++) hybridMap[hybridCells1[cell]] = cell; 1264 for (cell = 0; cell < hc2; cell++) hybridMap[hybridCells2[cell]] = cell + hc1; 1265 ierr = PetscFree(hybridCells1);CHKERRQ(ierr); 1266 ierr = PetscFree(hybridCells2);CHKERRQ(ierr); 1267 } 1268 1269 } 1270 1271 /* Allocate the cell-vertex mesh */ 1272 ierr = DMPlexSetChart(*dm, 0, trueNumCells+numVertices);CHKERRQ(ierr); 1273 for (cell = 0, c = 0; c < numCells; ++c) { 1274 if (gmsh_elem[c].dim == dim) { 1275 ierr = DMPlexSetConeSize(*dm, hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].numNodes);CHKERRQ(ierr); 1276 cell++; 1277 } 1278 } 1279 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1280 /* Add cell-vertex connections */ 1281 for (cell = 0, c = 0; c < numCells; ++c) { 1282 if (gmsh_elem[c].dim == dim) { 1283 PetscInt pcone[8], corner; 1284 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1285 const PetscInt cc = gmsh_elem[c].nodes[corner] - shift; 1286 pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + trueNumCells; 1287 } 1288 if (dim == 3) { 1289 /* Tetrahedra are inverted */ 1290 if (gmsh_elem[c].cellType == 4) { 1291 PetscInt tmp = pcone[0]; 1292 pcone[0] = pcone[1]; 1293 pcone[1] = tmp; 1294 } 1295 /* Hexahedra are inverted */ 1296 if (gmsh_elem[c].cellType == 5) { 1297 PetscInt tmp = pcone[1]; 1298 pcone[1] = pcone[3]; 1299 pcone[3] = tmp; 1300 } 1301 /* Prisms are inverted */ 1302 if (gmsh_elem[c].cellType == 6) { 1303 PetscInt tmp; 1304 1305 tmp = pcone[1]; 1306 pcone[1] = pcone[2]; 1307 pcone[2] = tmp; 1308 tmp = pcone[4]; 1309 pcone[4] = pcone[5]; 1310 pcone[5] = tmp; 1311 } 1312 } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */ 1313 /* quads are input to PLEX as prisms */ 1314 if (gmsh_elem[c].cellType == 3) { 1315 PetscInt tmp = pcone[2]; 1316 pcone[2] = pcone[3]; 1317 pcone[3] = tmp; 1318 } 1319 } 1320 ierr = DMPlexSetCone(*dm, hybridMap ? hybridMap[cell] : cell, pcone);CHKERRQ(ierr); 1321 cell++; 1322 } 1323 } 1324 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 1325 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1326 ierr = DMPlexSetHybridBounds(*dm, cMax, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1327 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1328 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1329 if (interpolate) { 1330 DM idm; 1331 1332 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1333 ierr = DMDestroy(dm);CHKERRQ(ierr); 1334 *dm = idm; 1335 } 1336 1337 if (usemarker && !interpolate && dim > 1) SETERRQ(comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation"); 1338 if (!rank && usemarker) { 1339 PetscInt f, fStart, fEnd; 1340 1341 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 1342 ierr = DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1343 for (f = fStart; f < fEnd; ++f) { 1344 PetscInt suppSize; 1345 1346 ierr = DMPlexGetSupportSize(*dm, f, &suppSize);CHKERRQ(ierr); 1347 if (suppSize == 1) { 1348 PetscInt *cone = NULL, coneSize, p; 1349 1350 ierr = DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);CHKERRQ(ierr); 1351 for (p = 0; p < coneSize; p += 2) { 1352 ierr = DMSetLabelValue(*dm, "marker", cone[p], 1);CHKERRQ(ierr); 1353 } 1354 ierr = DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);CHKERRQ(ierr); 1355 } 1356 } 1357 } 1358 1359 if (!rank) { 1360 PetscInt vStart, vEnd; 1361 1362 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1363 for (cell = 0, c = 0; c < numCells; ++c) { 1364 1365 /* Create face sets */ 1366 if (interpolate && gmsh_elem[c].dim == dim-1) { 1367 const PetscInt *join; 1368 PetscInt joinSize, pcone[8], corner; 1369 /* Find the relevant facet with vertex joins */ 1370 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1371 const PetscInt cc = gmsh_elem[c].nodes[corner] - shift; 1372 pcone[corner] = (periodicMap ? periodicMap[cc] : cc) + vStart; 1373 } 1374 ierr = DMPlexGetFullJoin(*dm, gmsh_elem[c].numNodes, pcone, &joinSize, &join);CHKERRQ(ierr); 1375 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); 1376 ierr = DMSetLabelValue(*dm, "Face Sets", join[0], gmsh_elem[c].tags[0]);CHKERRQ(ierr); 1377 ierr = DMPlexRestoreJoin(*dm, gmsh_elem[c].numNodes, (const PetscInt *) pcone, &joinSize, &join);CHKERRQ(ierr); 1378 } 1379 1380 /* Create cell sets */ 1381 if (gmsh_elem[c].dim == dim) { 1382 if (gmsh_elem[c].numTags > 0) { 1383 ierr = DMSetLabelValue(*dm, "Cell Sets", hybridMap ? hybridMap[cell] : cell, gmsh_elem[c].tags[0]);CHKERRQ(ierr); 1384 } 1385 cell++; 1386 } 1387 1388 /* Create vertex sets */ 1389 if (gmsh_elem[c].dim == 0) { 1390 if (gmsh_elem[c].numTags > 0) { 1391 const PetscInt cc = gmsh_elem[c].nodes[0] - shift; 1392 const PetscInt vid = (periodicMap ? periodicMap[cc] : cc) + vStart; 1393 ierr = DMSetLabelValue(*dm, "Vertex Sets", vid, gmsh_elem[c].tags[0]);CHKERRQ(ierr); 1394 } 1395 } 1396 } 1397 } 1398 1399 /* Create coordinates */ 1400 if (embedDim < 0) embedDim = dim; 1401 ierr = DMSetCoordinateDim(*dm, embedDim);CHKERRQ(ierr); 1402 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 1403 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1404 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr); 1405 if (periodicMap) { /* we need to localize coordinates on cells */ 1406 ierr = PetscSectionSetChart(coordSection, 0, trueNumCells + numVertices);CHKERRQ(ierr); 1407 } else { 1408 ierr = PetscSectionSetChart(coordSection, trueNumCells, trueNumCells + numVertices);CHKERRQ(ierr); 1409 } 1410 for (v = trueNumCells; v < trueNumCells+numVertices; ++v) { 1411 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr); 1412 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr); 1413 } 1414 if (periodicMap) { 1415 ierr = PetscBTCreate(trueNumCells, &periodicC);CHKERRQ(ierr); 1416 for (cell = 0, c = 0; c < numCells; ++c) { 1417 if (gmsh_elem[c].dim == dim) { 1418 PetscInt corner; 1419 PetscBool pc = PETSC_FALSE; 1420 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1421 pc = (PetscBool)(pc || PetscBTLookup(periodicV, gmsh_elem[c].nodes[corner] - shift)); 1422 } 1423 if (pc) { 1424 PetscInt dof = gmsh_elem[c].numNodes*embedDim; 1425 PetscInt ucell = hybridMap ? hybridMap[cell] : cell; 1426 ierr = PetscBTSet(periodicC, ucell);CHKERRQ(ierr); 1427 ierr = PetscSectionSetDof(coordSection, ucell, dof);CHKERRQ(ierr); 1428 ierr = PetscSectionSetFieldDof(coordSection, ucell, 0, dof);CHKERRQ(ierr); 1429 } 1430 cell++; 1431 } 1432 } 1433 } 1434 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1435 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 1436 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 1437 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1438 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 1439 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr); 1440 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 1441 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 1442 if (periodicMap) { 1443 PetscInt off; 1444 1445 for (cell = 0, c = 0; c < numCells; ++c) { 1446 PetscInt pcone[8], corner; 1447 if (gmsh_elem[c].dim == dim) { 1448 PetscInt ucell = hybridMap ? hybridMap[cell] : cell; 1449 if (PetscUnlikely(PetscBTLookup(periodicC, ucell))) { 1450 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1451 pcone[corner] = gmsh_elem[c].nodes[corner] - shift; 1452 } 1453 if (dim == 3) { 1454 /* Tetrahedra are inverted */ 1455 if (gmsh_elem[c].cellType == 4) { 1456 PetscInt tmp = pcone[0]; 1457 pcone[0] = pcone[1]; 1458 pcone[1] = tmp; 1459 } 1460 /* Hexahedra are inverted */ 1461 if (gmsh_elem[c].cellType == 5) { 1462 PetscInt tmp = pcone[1]; 1463 pcone[1] = pcone[3]; 1464 pcone[3] = tmp; 1465 } 1466 /* Prisms are inverted */ 1467 if (gmsh_elem[c].cellType == 6) { 1468 PetscInt tmp; 1469 1470 tmp = pcone[1]; 1471 pcone[1] = pcone[2]; 1472 pcone[2] = tmp; 1473 tmp = pcone[4]; 1474 pcone[4] = pcone[5]; 1475 pcone[5] = tmp; 1476 } 1477 } else if (dim == 2 && hybridMap && hybridMap[cell] >= cMax) { /* hybrid cells */ 1478 /* quads are input to PLEX as prisms */ 1479 if (gmsh_elem[c].cellType == 3) { 1480 PetscInt tmp = pcone[2]; 1481 pcone[2] = pcone[3]; 1482 pcone[3] = tmp; 1483 } 1484 } 1485 ierr = PetscSectionGetOffset(coordSection, ucell, &off);CHKERRQ(ierr); 1486 for (corner = 0; corner < gmsh_elem[c].numNodes; ++corner) { 1487 v = pcone[corner]; 1488 for (d = 0; d < embedDim; ++d) { 1489 coords[off++] = (PetscReal) coordsIn[v*3+d]; 1490 } 1491 } 1492 } 1493 cell++; 1494 } 1495 } 1496 for (v = 0; v < numVertices; ++v) { 1497 ierr = PetscSectionGetOffset(coordSection, v + trueNumCells, &off);CHKERRQ(ierr); 1498 for (d = 0; d < embedDim; ++d) { 1499 coords[off+d] = (PetscReal) coordsIn[periodicMapI[v]*3+d]; 1500 } 1501 } 1502 } else { 1503 for (v = 0; v < numVertices; ++v) { 1504 for (d = 0; d < embedDim; ++d) { 1505 coords[v*embedDim+d] = (PetscReal) coordsIn[v*3+d]; 1506 } 1507 } 1508 } 1509 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 1510 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1511 1512 periodic = periodicMap ? PETSC_TRUE : PETSC_FALSE; 1513 ierr = MPI_Bcast(&periodic, 1, MPIU_BOOL, 0, comm);CHKERRQ(ierr); 1514 ierr = DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);CHKERRQ(ierr); 1515 1516 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 1517 ierr = PetscFree(gmsh_elem);CHKERRQ(ierr); 1518 ierr = PetscFree(hybridMap);CHKERRQ(ierr); 1519 ierr = PetscFree(periodicMap);CHKERRQ(ierr); 1520 ierr = PetscFree(periodicMapI);CHKERRQ(ierr); 1521 ierr = PetscBTDestroy(&periodicV);CHKERRQ(ierr); 1522 ierr = PetscBTDestroy(&periodicC);CHKERRQ(ierr); 1523 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1524 1525 ierr = PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,0,0,0);CHKERRQ(ierr); 1526 PetscFunctionReturn(0); 1527 } 1528