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