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