1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petsc/private/hashmapi.h> 4 5 #include <../src/dm/impls/plex/gmshlex.h> 6 7 #define GMSH_LEXORDER_ITEM(T, p) \ 8 static int *GmshLexOrder_##T##_##p(void) \ 9 { \ 10 static int Gmsh_LexOrder_##T##_##p[GmshNumNodes_##T(p)] = {-1}; \ 11 int *lex = Gmsh_LexOrder_##T##_##p; \ 12 if (lex[0] == -1) (void)GmshLexOrder_##T(p, lex, 0); \ 13 return lex; \ 14 } 15 16 static int *GmshLexOrder_QUA_2_Serendipity(void) 17 { 18 static int Gmsh_LexOrder_QUA_2_Serendipity[9] = {-1}; 19 int *lex = Gmsh_LexOrder_QUA_2_Serendipity; 20 if (lex[0] == -1) { 21 /* Vertices */ 22 lex[0] = 0; lex[2] = 1; lex[8] = 2; lex[6] = 3; 23 /* Edges */ 24 lex[1] = 4; lex[5] = 5; lex[7] = 6; lex[3] = 7; 25 /* Cell */ 26 lex[4] = -8; 27 } 28 return lex; 29 } 30 31 static int *GmshLexOrder_HEX_2_Serendipity(void) 32 { 33 static int Gmsh_LexOrder_HEX_2_Serendipity[27] = {-1}; 34 int *lex = Gmsh_LexOrder_HEX_2_Serendipity; 35 if (lex[0] == -1) { 36 /* Vertices */ 37 lex[ 0] = 0; lex[ 2] = 1; lex[ 8] = 2; lex[ 6] = 3; 38 lex[18] = 4; lex[20] = 5; lex[26] = 6; lex[24] = 7; 39 /* Edges */ 40 lex[ 1] = 8; lex[ 3] = 9; lex[ 9] = 10; lex[ 5] = 11; 41 lex[11] = 12; lex[ 7] = 13; lex[17] = 14; lex[15] = 15; 42 lex[19] = 16; lex[21] = 17; lex[23] = 18; lex[25] = 19; 43 /* Faces */ 44 lex[ 4] = -20; lex[10] = -21; lex[12] = -22; lex[14] = -23; 45 lex[16] = -24; lex[22] = -25; 46 /* Cell */ 47 lex[13] = -26; 48 } 49 return lex; 50 } 51 52 #define GMSH_LEXORDER_LIST(T) \ 53 GMSH_LEXORDER_ITEM(T, 1) \ 54 GMSH_LEXORDER_ITEM(T, 2) \ 55 GMSH_LEXORDER_ITEM(T, 3) \ 56 GMSH_LEXORDER_ITEM(T, 4) \ 57 GMSH_LEXORDER_ITEM(T, 5) \ 58 GMSH_LEXORDER_ITEM(T, 6) \ 59 GMSH_LEXORDER_ITEM(T, 7) \ 60 GMSH_LEXORDER_ITEM(T, 8) \ 61 GMSH_LEXORDER_ITEM(T, 9) \ 62 GMSH_LEXORDER_ITEM(T, 10) 63 64 GMSH_LEXORDER_ITEM(VTX, 0) 65 GMSH_LEXORDER_LIST(SEG) 66 GMSH_LEXORDER_LIST(TRI) 67 GMSH_LEXORDER_LIST(QUA) 68 GMSH_LEXORDER_LIST(TET) 69 GMSH_LEXORDER_LIST(HEX) 70 GMSH_LEXORDER_LIST(PRI) 71 GMSH_LEXORDER_LIST(PYR) 72 73 typedef enum { 74 GMSH_VTX = 0, 75 GMSH_SEG = 1, 76 GMSH_TRI = 2, 77 GMSH_QUA = 3, 78 GMSH_TET = 4, 79 GMSH_HEX = 5, 80 GMSH_PRI = 6, 81 GMSH_PYR = 7, 82 GMSH_NUM_POLYTOPES = 8 83 } GmshPolytopeType; 84 85 typedef struct { 86 int cellType; 87 int polytope; 88 int dim; 89 int order; 90 int numVerts; 91 int numNodes; 92 int* (*lexorder)(void); 93 } GmshCellInfo; 94 95 #define GmshCellEntry(cellType, polytope, dim, order) \ 96 {cellType, GMSH_##polytope, dim, order, \ 97 GmshNumNodes_##polytope(1), \ 98 GmshNumNodes_##polytope(order), \ 99 GmshLexOrder_##polytope##_##order} 100 101 static const GmshCellInfo GmshCellTable[] = { 102 GmshCellEntry( 15, VTX, 0, 0), 103 104 GmshCellEntry( 1, SEG, 1, 1), 105 GmshCellEntry( 8, SEG, 1, 2), 106 GmshCellEntry( 26, SEG, 1, 3), 107 GmshCellEntry( 27, SEG, 1, 4), 108 GmshCellEntry( 28, SEG, 1, 5), 109 GmshCellEntry( 62, SEG, 1, 6), 110 GmshCellEntry( 63, SEG, 1, 7), 111 GmshCellEntry( 64, SEG, 1, 8), 112 GmshCellEntry( 65, SEG, 1, 9), 113 GmshCellEntry( 66, SEG, 1, 10), 114 115 GmshCellEntry( 2, TRI, 2, 1), 116 GmshCellEntry( 9, TRI, 2, 2), 117 GmshCellEntry( 21, TRI, 2, 3), 118 GmshCellEntry( 23, TRI, 2, 4), 119 GmshCellEntry( 25, TRI, 2, 5), 120 GmshCellEntry( 42, TRI, 2, 6), 121 GmshCellEntry( 43, TRI, 2, 7), 122 GmshCellEntry( 44, TRI, 2, 8), 123 GmshCellEntry( 45, TRI, 2, 9), 124 GmshCellEntry( 46, TRI, 2, 10), 125 126 GmshCellEntry( 3, QUA, 2, 1), 127 GmshCellEntry( 10, QUA, 2, 2), 128 {16, GMSH_QUA, 2, 2, 4, 8, GmshLexOrder_QUA_2_Serendipity}, 129 GmshCellEntry( 36, QUA, 2, 3), 130 GmshCellEntry( 37, QUA, 2, 4), 131 GmshCellEntry( 38, QUA, 2, 5), 132 GmshCellEntry( 47, QUA, 2, 6), 133 GmshCellEntry( 48, QUA, 2, 7), 134 GmshCellEntry( 49, QUA, 2, 8), 135 GmshCellEntry( 50, QUA, 2, 9), 136 GmshCellEntry( 51, QUA, 2, 10), 137 138 GmshCellEntry( 4, TET, 3, 1), 139 GmshCellEntry( 11, TET, 3, 2), 140 GmshCellEntry( 29, TET, 3, 3), 141 GmshCellEntry( 30, TET, 3, 4), 142 GmshCellEntry( 31, TET, 3, 5), 143 GmshCellEntry( 71, TET, 3, 6), 144 GmshCellEntry( 72, TET, 3, 7), 145 GmshCellEntry( 73, TET, 3, 8), 146 GmshCellEntry( 74, TET, 3, 9), 147 GmshCellEntry( 75, TET, 3, 10), 148 149 GmshCellEntry( 5, HEX, 3, 1), 150 GmshCellEntry( 12, HEX, 3, 2), 151 {17, GMSH_HEX, 3, 2, 8, 20, GmshLexOrder_HEX_2_Serendipity}, 152 GmshCellEntry( 92, HEX, 3, 3), 153 GmshCellEntry( 93, HEX, 3, 4), 154 GmshCellEntry( 94, HEX, 3, 5), 155 GmshCellEntry( 95, HEX, 3, 6), 156 GmshCellEntry( 96, HEX, 3, 7), 157 GmshCellEntry( 97, HEX, 3, 8), 158 GmshCellEntry( 98, HEX, 3, 9), 159 GmshCellEntry( -1, HEX, 3, 10), 160 161 GmshCellEntry( 6, PRI, 3, 1), 162 GmshCellEntry( 13, PRI, 3, 2), 163 GmshCellEntry( 90, PRI, 3, 3), 164 GmshCellEntry( 91, PRI, 3, 4), 165 GmshCellEntry(106, PRI, 3, 5), 166 GmshCellEntry(107, PRI, 3, 6), 167 GmshCellEntry(108, PRI, 3, 7), 168 GmshCellEntry(109, PRI, 3, 8), 169 GmshCellEntry(110, PRI, 3, 9), 170 GmshCellEntry( -1, PRI, 3, 10), 171 172 GmshCellEntry( 7, PYR, 3, 1), 173 GmshCellEntry( 14, PYR, 3, 2), 174 GmshCellEntry(118, PYR, 3, 3), 175 GmshCellEntry(119, PYR, 3, 4), 176 GmshCellEntry(120, PYR, 3, 5), 177 GmshCellEntry(121, PYR, 3, 6), 178 GmshCellEntry(122, PYR, 3, 7), 179 GmshCellEntry(123, PYR, 3, 8), 180 GmshCellEntry(124, PYR, 3, 9), 181 GmshCellEntry( -1, PYR, 3, 10) 182 183 #if 0 184 {20, GMSH_TRI, 2, 3, 3, 9, NULL}, 185 {18, GMSH_PRI, 3, 2, 6, 15, NULL}, 186 {19, GMSH_PYR, 3, 2, 5, 13, NULL}, 187 #endif 188 }; 189 190 static GmshCellInfo GmshCellMap[150]; 191 192 static PetscErrorCode GmshCellInfoSetUp(void) 193 { 194 size_t i, n; 195 static PetscBool called = PETSC_FALSE; 196 197 if (called) return 0; 198 PetscFunctionBegin; 199 called = PETSC_TRUE; 200 n = sizeof(GmshCellMap)/sizeof(GmshCellMap[0]); 201 for (i = 0; i < n; ++i) { 202 GmshCellMap[i].cellType = -1; 203 GmshCellMap[i].polytope = -1; 204 } 205 n = sizeof(GmshCellTable)/sizeof(GmshCellTable[0]); 206 for (i = 0; i < n; ++i) { 207 if (GmshCellTable[i].cellType <= 0) continue; 208 GmshCellMap[GmshCellTable[i].cellType] = GmshCellTable[i]; 209 } 210 PetscFunctionReturn(0); 211 } 212 213 #define GmshCellTypeCheck(ct) 0; do { \ 214 const int _ct_ = (int)ct; \ 215 if (_ct_ < 0 || _ct_ >= (int)(sizeof(GmshCellMap)/sizeof(GmshCellMap[0]))) \ 216 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid Gmsh element type %d", _ct_); \ 217 if (GmshCellMap[_ct_].cellType != _ct_) \ 218 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported Gmsh element type %d", _ct_); \ 219 if (GmshCellMap[_ct_].polytope == -1) \ 220 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported Gmsh element type %d", _ct_); \ 221 } while (0) 222 223 typedef struct { 224 PetscViewer viewer; 225 int fileFormat; 226 int dataSize; 227 PetscBool binary; 228 PetscBool byteSwap; 229 size_t wlen; 230 void *wbuf; 231 size_t slen; 232 void *sbuf; 233 PetscInt *nbuf; 234 PetscInt nodeStart; 235 PetscInt nodeEnd; 236 PetscInt *nodeMap; 237 } GmshFile; 238 239 static PetscErrorCode GmshBufferGet(GmshFile *gmsh, size_t count, size_t eltsize, void *buf) 240 { 241 size_t size = count * eltsize; 242 PetscErrorCode ierr; 243 244 PetscFunctionBegin; 245 if (gmsh->wlen < size) { 246 ierr = PetscFree(gmsh->wbuf);CHKERRQ(ierr); 247 ierr = PetscMalloc(size, &gmsh->wbuf);CHKERRQ(ierr); 248 gmsh->wlen = size; 249 } 250 *(void**)buf = size ? gmsh->wbuf : NULL; 251 PetscFunctionReturn(0); 252 } 253 254 static PetscErrorCode GmshBufferSizeGet(GmshFile *gmsh, size_t count, void *buf) 255 { 256 size_t dataSize = (size_t)gmsh->dataSize; 257 size_t size = count * dataSize; 258 PetscErrorCode ierr; 259 260 PetscFunctionBegin; 261 if (gmsh->slen < size) { 262 ierr = PetscFree(gmsh->sbuf);CHKERRQ(ierr); 263 ierr = PetscMalloc(size, &gmsh->sbuf);CHKERRQ(ierr); 264 gmsh->slen = size; 265 } 266 *(void**)buf = size ? gmsh->sbuf : NULL; 267 PetscFunctionReturn(0); 268 } 269 270 static PetscErrorCode GmshRead(GmshFile *gmsh, void *buf, PetscInt count, PetscDataType dtype) 271 { 272 PetscErrorCode ierr; 273 PetscFunctionBegin; 274 ierr = PetscViewerRead(gmsh->viewer, buf, count, NULL, dtype);CHKERRQ(ierr); 275 if (gmsh->byteSwap) {ierr = PetscByteSwap(buf, dtype, count);CHKERRQ(ierr);} 276 PetscFunctionReturn(0); 277 } 278 279 static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscInt count) 280 { 281 PetscErrorCode ierr; 282 PetscFunctionBegin; 283 ierr = PetscViewerRead(gmsh->viewer, buf, count, NULL, PETSC_STRING);CHKERRQ(ierr); 284 PetscFunctionReturn(0); 285 } 286 287 static PetscErrorCode GmshMatch(PETSC_UNUSED GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN], PetscBool *match) 288 { 289 PetscErrorCode ierr; 290 PetscFunctionBegin; 291 ierr = PetscStrcmp(line, Section, match);CHKERRQ(ierr); 292 PetscFunctionReturn(0); 293 } 294 295 static PetscErrorCode GmshExpect(GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN]) 296 { 297 PetscBool match; 298 PetscErrorCode ierr; 299 300 PetscFunctionBegin; 301 ierr = GmshMatch(gmsh, Section, line, &match);CHKERRQ(ierr); 302 PetscCheckFalse(!match,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file, expecting %s",Section); 303 PetscFunctionReturn(0); 304 } 305 306 static PetscErrorCode GmshReadSection(GmshFile *gmsh, char line[PETSC_MAX_PATH_LEN]) 307 { 308 PetscBool match; 309 PetscErrorCode ierr; 310 311 PetscFunctionBegin; 312 while (PETSC_TRUE) { 313 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 314 ierr = GmshMatch(gmsh, "$Comments", line, &match);CHKERRQ(ierr); 315 if (!match) break; 316 while (PETSC_TRUE) { 317 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 318 ierr = GmshMatch(gmsh, "$EndComments", line, &match);CHKERRQ(ierr); 319 if (match) break; 320 } 321 } 322 PetscFunctionReturn(0); 323 } 324 325 static PetscErrorCode GmshReadEndSection(GmshFile *gmsh, const char EndSection[], char line[PETSC_MAX_PATH_LEN]) 326 { 327 PetscErrorCode ierr; 328 PetscFunctionBegin; 329 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 330 ierr = GmshExpect(gmsh, EndSection, line);CHKERRQ(ierr); 331 PetscFunctionReturn(0); 332 } 333 334 static PetscErrorCode GmshReadSize(GmshFile *gmsh, PetscInt *buf, PetscInt count) 335 { 336 PetscInt i; 337 size_t dataSize = (size_t)gmsh->dataSize; 338 PetscErrorCode ierr; 339 340 PetscFunctionBegin; 341 if (dataSize == sizeof(PetscInt)) { 342 ierr = GmshRead(gmsh, buf, count, PETSC_INT);CHKERRQ(ierr); 343 } else if (dataSize == sizeof(int)) { 344 int *ibuf = NULL; 345 ierr = GmshBufferSizeGet(gmsh, count, &ibuf);CHKERRQ(ierr); 346 ierr = GmshRead(gmsh, ibuf, count, PETSC_ENUM);CHKERRQ(ierr); 347 for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i]; 348 } else if (dataSize == sizeof(long)) { 349 long *ibuf = NULL; 350 ierr = GmshBufferSizeGet(gmsh, count, &ibuf);CHKERRQ(ierr); 351 ierr = GmshRead(gmsh, ibuf, count, PETSC_LONG);CHKERRQ(ierr); 352 for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i]; 353 } else if (dataSize == sizeof(PetscInt64)) { 354 PetscInt64 *ibuf = NULL; 355 ierr = GmshBufferSizeGet(gmsh, count, &ibuf);CHKERRQ(ierr); 356 ierr = GmshRead(gmsh, ibuf, count, PETSC_INT64);CHKERRQ(ierr); 357 for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i]; 358 } 359 PetscFunctionReturn(0); 360 } 361 362 static PetscErrorCode GmshReadInt(GmshFile *gmsh, int *buf, PetscInt count) 363 { 364 PetscErrorCode ierr; 365 PetscFunctionBegin; 366 ierr = GmshRead(gmsh, buf, count, PETSC_ENUM);CHKERRQ(ierr); 367 PetscFunctionReturn(0); 368 } 369 370 static PetscErrorCode GmshReadDouble(GmshFile *gmsh, double *buf, PetscInt count) 371 { 372 PetscErrorCode ierr; 373 PetscFunctionBegin; 374 ierr = GmshRead(gmsh, buf, count, PETSC_DOUBLE);CHKERRQ(ierr); 375 PetscFunctionReturn(0); 376 } 377 378 typedef struct { 379 PetscInt id; /* Entity ID */ 380 PetscInt dim; /* Dimension */ 381 double bbox[6]; /* Bounding box */ 382 PetscInt numTags; /* Size of tag array */ 383 int tags[4]; /* Tag array */ 384 } GmshEntity; 385 386 typedef struct { 387 GmshEntity *entity[4]; 388 PetscHMapI entityMap[4]; 389 } GmshEntities; 390 391 static PetscErrorCode GmshEntitiesCreate(PetscInt count[4], GmshEntities **entities) 392 { 393 PetscInt dim; 394 PetscErrorCode ierr; 395 396 PetscFunctionBegin; 397 ierr = PetscNew(entities);CHKERRQ(ierr); 398 for (dim = 0; dim < 4; ++dim) { 399 ierr = PetscCalloc1(count[dim], &(*entities)->entity[dim]);CHKERRQ(ierr); 400 ierr = PetscHMapICreate(&(*entities)->entityMap[dim]);CHKERRQ(ierr); 401 } 402 PetscFunctionReturn(0); 403 } 404 405 static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities) 406 { 407 PetscInt dim; 408 PetscErrorCode ierr; 409 410 PetscFunctionBegin; 411 if (!*entities) PetscFunctionReturn(0); 412 for (dim = 0; dim < 4; ++dim) { 413 ierr = PetscFree((*entities)->entity[dim]);CHKERRQ(ierr); 414 ierr = PetscHMapIDestroy(&(*entities)->entityMap[dim]);CHKERRQ(ierr); 415 } 416 ierr = PetscFree((*entities));CHKERRQ(ierr); 417 PetscFunctionReturn(0); 418 } 419 420 static PetscErrorCode GmshEntitiesAdd(GmshEntities *entities, PetscInt index, PetscInt dim, PetscInt eid, GmshEntity** entity) 421 { 422 PetscErrorCode ierr; 423 424 PetscFunctionBegin; 425 ierr = PetscHMapISet(entities->entityMap[dim], eid, index);CHKERRQ(ierr); 426 entities->entity[dim][index].dim = dim; 427 entities->entity[dim][index].id = eid; 428 if (entity) *entity = &entities->entity[dim][index]; 429 PetscFunctionReturn(0); 430 } 431 432 static PetscErrorCode GmshEntitiesGet(GmshEntities *entities, PetscInt dim, PetscInt eid, GmshEntity** entity) 433 { 434 PetscInt index; 435 PetscErrorCode ierr; 436 437 PetscFunctionBegin; 438 ierr = PetscHMapIGet(entities->entityMap[dim], eid, &index);CHKERRQ(ierr); 439 *entity = &entities->entity[dim][index]; 440 PetscFunctionReturn(0); 441 } 442 443 typedef struct { 444 PetscInt *id; /* Node IDs */ 445 double *xyz; /* Coordinates */ 446 } GmshNodes; 447 448 static PetscErrorCode GmshNodesCreate(PetscInt count, GmshNodes **nodes) 449 { 450 PetscErrorCode ierr; 451 452 PetscFunctionBegin; 453 ierr = PetscNew(nodes);CHKERRQ(ierr); 454 ierr = PetscMalloc1(count*1, &(*nodes)->id);CHKERRQ(ierr); 455 ierr = PetscMalloc1(count*3, &(*nodes)->xyz);CHKERRQ(ierr); 456 PetscFunctionReturn(0); 457 } 458 459 static PetscErrorCode GmshNodesDestroy(GmshNodes **nodes) 460 { 461 PetscErrorCode ierr; 462 PetscFunctionBegin; 463 if (!*nodes) PetscFunctionReturn(0); 464 ierr = PetscFree((*nodes)->id);CHKERRQ(ierr); 465 ierr = PetscFree((*nodes)->xyz);CHKERRQ(ierr); 466 ierr = PetscFree((*nodes));CHKERRQ(ierr); 467 PetscFunctionReturn(0); 468 } 469 470 typedef struct { 471 PetscInt id; /* Element ID */ 472 PetscInt dim; /* Dimension */ 473 PetscInt cellType; /* Cell type */ 474 PetscInt numVerts; /* Size of vertex array */ 475 PetscInt numNodes; /* Size of node array */ 476 PetscInt *nodes; /* Vertex/Node array */ 477 PetscInt numTags; /* Size of tag array */ 478 int tags[4]; /* Tag array */ 479 } GmshElement; 480 481 static PetscErrorCode GmshElementsCreate(PetscInt count, GmshElement **elements) 482 { 483 PetscErrorCode ierr; 484 485 PetscFunctionBegin; 486 ierr = PetscCalloc1(count, elements);CHKERRQ(ierr); 487 PetscFunctionReturn(0); 488 } 489 490 static PetscErrorCode GmshElementsDestroy(GmshElement **elements) 491 { 492 PetscErrorCode ierr; 493 494 PetscFunctionBegin; 495 if (!*elements) PetscFunctionReturn(0); 496 ierr = PetscFree(*elements);CHKERRQ(ierr); 497 PetscFunctionReturn(0); 498 } 499 500 typedef struct { 501 PetscInt dim; 502 PetscInt order; 503 GmshEntities *entities; 504 PetscInt numNodes; 505 GmshNodes *nodelist; 506 PetscInt numElems; 507 GmshElement *elements; 508 PetscInt numVerts; 509 PetscInt numCells; 510 PetscInt *periodMap; 511 PetscInt *vertexMap; 512 PetscSegBuffer segbuf; 513 PetscInt numRegions; 514 PetscInt *regionTags; 515 char **regionNames; 516 } GmshMesh; 517 518 static PetscErrorCode GmshMeshCreate(GmshMesh **mesh) 519 { 520 PetscErrorCode ierr; 521 522 PetscFunctionBegin; 523 ierr = PetscNew(mesh);CHKERRQ(ierr); 524 ierr = PetscSegBufferCreate(sizeof(PetscInt), 0, &(*mesh)->segbuf);CHKERRQ(ierr); 525 PetscFunctionReturn(0); 526 } 527 528 static PetscErrorCode GmshMeshDestroy(GmshMesh **mesh) 529 { 530 PetscInt r; 531 PetscErrorCode ierr; 532 533 PetscFunctionBegin; 534 if (!*mesh) PetscFunctionReturn(0); 535 ierr = GmshEntitiesDestroy(&(*mesh)->entities);CHKERRQ(ierr); 536 ierr = GmshNodesDestroy(&(*mesh)->nodelist);CHKERRQ(ierr); 537 ierr = GmshElementsDestroy(&(*mesh)->elements);CHKERRQ(ierr); 538 ierr = PetscFree((*mesh)->periodMap);CHKERRQ(ierr); 539 ierr = PetscFree((*mesh)->vertexMap);CHKERRQ(ierr); 540 ierr = PetscSegBufferDestroy(&(*mesh)->segbuf);CHKERRQ(ierr); 541 for (r = 0; r < (*mesh)->numRegions; ++r) {ierr = PetscFree((*mesh)->regionNames[r]);CHKERRQ(ierr);} 542 ierr = PetscFree2((*mesh)->regionTags, (*mesh)->regionNames);CHKERRQ(ierr); 543 ierr = PetscFree((*mesh));CHKERRQ(ierr); 544 PetscFunctionReturn(0); 545 } 546 547 static PetscErrorCode GmshReadNodes_v22(GmshFile *gmsh, GmshMesh *mesh) 548 { 549 PetscViewer viewer = gmsh->viewer; 550 PetscBool byteSwap = gmsh->byteSwap; 551 char line[PETSC_MAX_PATH_LEN]; 552 int n, num, nid, snum; 553 GmshNodes *nodes; 554 PetscErrorCode ierr; 555 556 PetscFunctionBegin; 557 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 558 snum = sscanf(line, "%d", &num); 559 PetscCheckFalse(snum != 1,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 560 ierr = GmshNodesCreate(num, &nodes);CHKERRQ(ierr); 561 mesh->numNodes = num; 562 mesh->nodelist = nodes; 563 for (n = 0; n < num; ++n) { 564 double *xyz = nodes->xyz + n*3; 565 ierr = PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 566 ierr = PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);CHKERRQ(ierr); 567 if (byteSwap) {ierr = PetscByteSwap(&nid, PETSC_ENUM, 1);CHKERRQ(ierr);} 568 if (byteSwap) {ierr = PetscByteSwap(xyz, PETSC_DOUBLE, 3);CHKERRQ(ierr);} 569 nodes->id[n] = nid; 570 } 571 PetscFunctionReturn(0); 572 } 573 574 /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the 575 file contents multiple times to figure out the true number of cells and facets 576 in the given mesh. To make this more efficient we read the file contents only 577 once and store them in memory, while determining the true number of cells. */ 578 static PetscErrorCode GmshReadElements_v22(GmshFile* gmsh, GmshMesh *mesh) 579 { 580 PetscViewer viewer = gmsh->viewer; 581 PetscBool binary = gmsh->binary; 582 PetscBool byteSwap = gmsh->byteSwap; 583 char line[PETSC_MAX_PATH_LEN]; 584 int i, c, p, num, ibuf[1+4+1000], snum; 585 int cellType, numElem, numVerts, numNodes, numTags; 586 GmshElement *elements; 587 PetscInt *nodeMap = gmsh->nodeMap; 588 PetscErrorCode ierr; 589 590 PetscFunctionBegin; 591 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 592 snum = sscanf(line, "%d", &num); 593 PetscCheckFalse(snum != 1,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 594 ierr = GmshElementsCreate(num, &elements);CHKERRQ(ierr); 595 mesh->numElems = num; 596 mesh->elements = elements; 597 for (c = 0; c < num;) { 598 ierr = PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);CHKERRQ(ierr); 599 if (byteSwap) {ierr = PetscByteSwap(ibuf, PETSC_ENUM, 3);CHKERRQ(ierr);} 600 601 cellType = binary ? ibuf[0] : ibuf[1]; 602 numElem = binary ? ibuf[1] : 1; 603 numTags = ibuf[2]; 604 605 ierr = GmshCellTypeCheck(cellType);CHKERRQ(ierr); 606 numVerts = GmshCellMap[cellType].numVerts; 607 numNodes = GmshCellMap[cellType].numNodes; 608 609 for (i = 0; i < numElem; ++i, ++c) { 610 GmshElement *element = elements + c; 611 const int off = binary ? 0 : 1, nint = 1 + numTags + numNodes - off; 612 const int *id = ibuf, *nodes = ibuf + 1 + numTags, *tags = ibuf + 1; 613 ierr = PetscViewerRead(viewer, ibuf+off, nint, NULL, PETSC_ENUM);CHKERRQ(ierr); 614 if (byteSwap) {ierr = PetscByteSwap(ibuf+off, PETSC_ENUM, nint);CHKERRQ(ierr);} 615 element->id = id[0]; 616 element->dim = GmshCellMap[cellType].dim; 617 element->cellType = cellType; 618 element->numVerts = numVerts; 619 element->numNodes = numNodes; 620 element->numTags = PetscMin(numTags, 4); 621 ierr = PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);CHKERRQ(ierr); 622 for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]]; 623 for (p = 0; p < element->numTags; p++) element->tags[p] = tags[p]; 624 } 625 } 626 PetscFunctionReturn(0); 627 } 628 629 /* 630 $Entities 631 numPoints(unsigned long) numCurves(unsigned long) 632 numSurfaces(unsigned long) numVolumes(unsigned long) 633 // points 634 tag(int) boxMinX(double) boxMinY(double) boxMinZ(double) 635 boxMaxX(double) boxMaxY(double) boxMaxZ(double) 636 numPhysicals(unsigned long) phyisicalTag[...](int) 637 ... 638 // curves 639 tag(int) boxMinX(double) boxMinY(double) boxMinZ(double) 640 boxMaxX(double) boxMaxY(double) boxMaxZ(double) 641 numPhysicals(unsigned long) physicalTag[...](int) 642 numBREPVert(unsigned long) tagBREPVert[...](int) 643 ... 644 // surfaces 645 tag(int) boxMinX(double) boxMinY(double) boxMinZ(double) 646 boxMaxX(double) boxMaxY(double) boxMaxZ(double) 647 numPhysicals(unsigned long) physicalTag[...](int) 648 numBREPCurve(unsigned long) tagBREPCurve[...](int) 649 ... 650 // volumes 651 tag(int) boxMinX(double) boxMinY(double) boxMinZ(double) 652 boxMaxX(double) boxMaxY(double) boxMaxZ(double) 653 numPhysicals(unsigned long) physicalTag[...](int) 654 numBREPSurfaces(unsigned long) tagBREPSurfaces[...](int) 655 ... 656 $EndEntities 657 */ 658 static PetscErrorCode GmshReadEntities_v40(GmshFile *gmsh, GmshMesh *mesh) 659 { 660 PetscViewer viewer = gmsh->viewer; 661 PetscBool byteSwap = gmsh->byteSwap; 662 long index, num, lbuf[4]; 663 int dim, eid, numTags, *ibuf, t; 664 PetscInt count[4], i; 665 GmshEntity *entity = NULL; 666 PetscErrorCode ierr; 667 668 PetscFunctionBegin; 669 ierr = PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG);CHKERRQ(ierr); 670 if (byteSwap) {ierr = PetscByteSwap(lbuf, PETSC_LONG, 4);CHKERRQ(ierr);} 671 for (i = 0; i < 4; ++i) count[i] = lbuf[i]; 672 ierr = GmshEntitiesCreate(count, &mesh->entities);CHKERRQ(ierr); 673 for (dim = 0; dim < 4; ++dim) { 674 for (index = 0; index < count[dim]; ++index) { 675 ierr = PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 676 if (byteSwap) {ierr = PetscByteSwap(&eid, PETSC_ENUM, 1);CHKERRQ(ierr);} 677 ierr = GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity);CHKERRQ(ierr); 678 ierr = PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE);CHKERRQ(ierr); 679 if (byteSwap) {ierr = PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6);CHKERRQ(ierr);} 680 ierr = PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 681 if (byteSwap) {ierr = PetscByteSwap(&num, PETSC_LONG, 1);CHKERRQ(ierr);} 682 ierr = GmshBufferGet(gmsh, num, sizeof(int), &ibuf);CHKERRQ(ierr); 683 ierr = PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);CHKERRQ(ierr); 684 if (byteSwap) {ierr = PetscByteSwap(ibuf, PETSC_ENUM, num);CHKERRQ(ierr);} 685 entity->numTags = numTags = (int) PetscMin(num, 4); 686 for (t = 0; t < numTags; ++t) entity->tags[t] = ibuf[t]; 687 if (dim == 0) continue; 688 ierr = PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 689 if (byteSwap) {ierr = PetscByteSwap(&num, PETSC_LONG, 1);CHKERRQ(ierr);} 690 ierr = GmshBufferGet(gmsh, num, sizeof(int), &ibuf);CHKERRQ(ierr); 691 ierr = PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);CHKERRQ(ierr); 692 if (byteSwap) {ierr = PetscByteSwap(ibuf, PETSC_ENUM, num);CHKERRQ(ierr);} 693 } 694 } 695 PetscFunctionReturn(0); 696 } 697 698 /* 699 $Nodes 700 numEntityBlocks(unsigned long) numNodes(unsigned long) 701 tagEntity(int) dimEntity(int) typeNode(int) numNodes(unsigned long) 702 tag(int) x(double) y(double) z(double) 703 ... 704 ... 705 $EndNodes 706 */ 707 static PetscErrorCode GmshReadNodes_v40(GmshFile *gmsh, GmshMesh *mesh) 708 { 709 PetscViewer viewer = gmsh->viewer; 710 PetscBool byteSwap = gmsh->byteSwap; 711 long block, node, n, numEntityBlocks, numTotalNodes, numNodes; 712 int info[3], nid; 713 GmshNodes *nodes; 714 PetscErrorCode ierr; 715 716 PetscFunctionBegin; 717 ierr = PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 718 if (byteSwap) {ierr = PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);CHKERRQ(ierr);} 719 ierr = PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 720 if (byteSwap) {ierr = PetscByteSwap(&numTotalNodes, PETSC_LONG, 1);CHKERRQ(ierr);} 721 ierr = GmshNodesCreate(numTotalNodes, &nodes);CHKERRQ(ierr); 722 mesh->numNodes = numTotalNodes; 723 mesh->nodelist = nodes; 724 for (n = 0, block = 0; block < numEntityBlocks; ++block) { 725 ierr = PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);CHKERRQ(ierr); 726 ierr = PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 727 if (byteSwap) {ierr = PetscByteSwap(&numNodes, PETSC_LONG, 1);CHKERRQ(ierr);} 728 if (gmsh->binary) { 729 size_t nbytes = sizeof(int) + 3*sizeof(double); 730 char *cbuf = NULL; /* dummy value to prevent warning from compiler about possible unitilized value */ 731 ierr = GmshBufferGet(gmsh, numNodes, nbytes, &cbuf);CHKERRQ(ierr); 732 ierr = PetscViewerRead(viewer, cbuf, numNodes*nbytes, NULL, PETSC_CHAR);CHKERRQ(ierr); 733 for (node = 0; node < numNodes; ++node, ++n) { 734 char *cnid = cbuf + node*nbytes, *cxyz = cnid + sizeof(int); 735 double *xyz = nodes->xyz + n*3; 736 if (!PetscBinaryBigEndian()) {ierr = PetscByteSwap(cnid, PETSC_ENUM, 1);CHKERRQ(ierr);} 737 if (!PetscBinaryBigEndian()) {ierr = PetscByteSwap(cxyz, PETSC_DOUBLE, 3);CHKERRQ(ierr);} 738 ierr = PetscMemcpy(&nid, cnid, sizeof(int));CHKERRQ(ierr); 739 ierr = PetscMemcpy(xyz, cxyz, 3*sizeof(double));CHKERRQ(ierr); 740 if (byteSwap) {ierr = PetscByteSwap(&nid, PETSC_ENUM, 1);CHKERRQ(ierr);} 741 if (byteSwap) {ierr = PetscByteSwap(xyz, PETSC_DOUBLE, 3);CHKERRQ(ierr);} 742 nodes->id[n] = nid; 743 } 744 } else { 745 for (node = 0; node < numNodes; ++node, ++n) { 746 double *xyz = nodes->xyz + n*3; 747 ierr = PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 748 ierr = PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);CHKERRQ(ierr); 749 if (byteSwap) {ierr = PetscByteSwap(&nid, PETSC_ENUM, 1);CHKERRQ(ierr);} 750 if (byteSwap) {ierr = PetscByteSwap(xyz, PETSC_DOUBLE, 3);CHKERRQ(ierr);} 751 nodes->id[n] = nid; 752 } 753 } 754 } 755 PetscFunctionReturn(0); 756 } 757 758 /* 759 $Elements 760 numEntityBlocks(unsigned long) numElements(unsigned long) 761 tagEntity(int) dimEntity(int) typeEle(int) numElements(unsigned long) 762 tag(int) numVert[...](int) 763 ... 764 ... 765 $EndElements 766 */ 767 static PetscErrorCode GmshReadElements_v40(GmshFile *gmsh, GmshMesh *mesh) 768 { 769 PetscViewer viewer = gmsh->viewer; 770 PetscBool byteSwap = gmsh->byteSwap; 771 long c, block, numEntityBlocks, numTotalElements, elem, numElements; 772 int p, info[3], *ibuf = NULL; 773 int eid, dim, cellType, numVerts, numNodes, numTags; 774 GmshEntity *entity = NULL; 775 GmshElement *elements; 776 PetscInt *nodeMap = gmsh->nodeMap; 777 PetscErrorCode ierr; 778 779 PetscFunctionBegin; 780 ierr = PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 781 if (byteSwap) {ierr = PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);CHKERRQ(ierr);} 782 ierr = PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 783 if (byteSwap) {ierr = PetscByteSwap(&numTotalElements, PETSC_LONG, 1);CHKERRQ(ierr);} 784 ierr = GmshElementsCreate(numTotalElements, &elements);CHKERRQ(ierr); 785 mesh->numElems = numTotalElements; 786 mesh->elements = elements; 787 for (c = 0, block = 0; block < numEntityBlocks; ++block) { 788 ierr = PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);CHKERRQ(ierr); 789 if (byteSwap) {ierr = PetscByteSwap(info, PETSC_ENUM, 3);CHKERRQ(ierr);} 790 eid = info[0]; dim = info[1]; cellType = info[2]; 791 ierr = GmshEntitiesGet(mesh->entities, dim, eid, &entity);CHKERRQ(ierr); 792 ierr = GmshCellTypeCheck(cellType);CHKERRQ(ierr); 793 numVerts = GmshCellMap[cellType].numVerts; 794 numNodes = GmshCellMap[cellType].numNodes; 795 numTags = entity->numTags; 796 ierr = PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 797 if (byteSwap) {ierr = PetscByteSwap(&numElements, PETSC_LONG, 1);CHKERRQ(ierr);} 798 ierr = GmshBufferGet(gmsh, (1+numNodes)*numElements, sizeof(int), &ibuf);CHKERRQ(ierr); 799 ierr = PetscViewerRead(viewer, ibuf, (1+numNodes)*numElements, NULL, PETSC_ENUM);CHKERRQ(ierr); 800 if (byteSwap) {ierr = PetscByteSwap(ibuf, PETSC_ENUM, (1+numNodes)*numElements);CHKERRQ(ierr);} 801 for (elem = 0; elem < numElements; ++elem, ++c) { 802 GmshElement *element = elements + c; 803 const int *id = ibuf + elem*(1+numNodes), *nodes = id + 1; 804 element->id = id[0]; 805 element->dim = dim; 806 element->cellType = cellType; 807 element->numVerts = numVerts; 808 element->numNodes = numNodes; 809 element->numTags = numTags; 810 ierr = PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);CHKERRQ(ierr); 811 for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]]; 812 for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p]; 813 } 814 } 815 PetscFunctionReturn(0); 816 } 817 818 static PetscErrorCode GmshReadPeriodic_v40(GmshFile *gmsh, PetscInt periodicMap[]) 819 { 820 PetscViewer viewer = gmsh->viewer; 821 int fileFormat = gmsh->fileFormat; 822 PetscBool binary = gmsh->binary; 823 PetscBool byteSwap = gmsh->byteSwap; 824 int numPeriodic, snum, i; 825 char line[PETSC_MAX_PATH_LEN]; 826 PetscInt *nodeMap = gmsh->nodeMap; 827 PetscErrorCode ierr; 828 829 PetscFunctionBegin; 830 if (fileFormat == 22 || !binary) { 831 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 832 snum = sscanf(line, "%d", &numPeriodic); 833 PetscCheckFalse(snum != 1,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 834 } else { 835 ierr = PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM);CHKERRQ(ierr); 836 if (byteSwap) {ierr = PetscByteSwap(&numPeriodic, PETSC_ENUM, 1);CHKERRQ(ierr);} 837 } 838 for (i = 0; i < numPeriodic; i++) { 839 int ibuf[3], correspondingDim = -1, correspondingTag = -1, primaryTag = -1, correspondingNode, primaryNode; 840 long j, nNodes; 841 double affine[16]; 842 843 if (fileFormat == 22 || !binary) { 844 ierr = PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);CHKERRQ(ierr); 845 snum = sscanf(line, "%d %d %d", &correspondingDim, &correspondingTag, &primaryTag); 846 PetscCheckFalse(snum != 3,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 847 } else { 848 ierr = PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);CHKERRQ(ierr); 849 if (byteSwap) {ierr = PetscByteSwap(ibuf, PETSC_ENUM, 3);CHKERRQ(ierr);} 850 correspondingDim = ibuf[0]; correspondingTag = ibuf[1]; primaryTag = ibuf[2]; 851 } 852 (void)correspondingDim; (void)correspondingTag; (void)primaryTag; /* unused */ 853 854 if (fileFormat == 22 || !binary) { 855 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 856 snum = sscanf(line, "%ld", &nNodes); 857 if (snum != 1) { /* discard transformation and try again */ 858 ierr = PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);CHKERRQ(ierr); 859 ierr = PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);CHKERRQ(ierr); 860 snum = sscanf(line, "%ld", &nNodes); 861 PetscCheckFalse(snum != 1,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 862 } 863 } else { 864 ierr = PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 865 if (byteSwap) {ierr = PetscByteSwap(&nNodes, PETSC_LONG, 1);CHKERRQ(ierr);} 866 if (nNodes == -1) { /* discard transformation and try again */ 867 ierr = PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE);CHKERRQ(ierr); 868 ierr = PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);CHKERRQ(ierr); 869 if (byteSwap) {ierr = PetscByteSwap(&nNodes, PETSC_LONG, 1);CHKERRQ(ierr);} 870 } 871 } 872 873 for (j = 0; j < nNodes; j++) { 874 if (fileFormat == 22 || !binary) { 875 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr); 876 snum = sscanf(line, "%d %d", &correspondingNode, &primaryNode); 877 PetscCheckFalse(snum != 2,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 878 } else { 879 ierr = PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM);CHKERRQ(ierr); 880 if (byteSwap) {ierr = PetscByteSwap(ibuf, PETSC_ENUM, 2);CHKERRQ(ierr);} 881 correspondingNode = ibuf[0]; primaryNode = ibuf[1]; 882 } 883 correspondingNode = (int) nodeMap[correspondingNode]; 884 primaryNode = (int) nodeMap[primaryNode]; 885 periodicMap[correspondingNode] = primaryNode; 886 } 887 } 888 PetscFunctionReturn(0); 889 } 890 891 /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format 892 $Entities 893 numPoints(size_t) numCurves(size_t) 894 numSurfaces(size_t) numVolumes(size_t) 895 pointTag(int) X(double) Y(double) Z(double) 896 numPhysicalTags(size_t) physicalTag(int) ... 897 ... 898 curveTag(int) minX(double) minY(double) minZ(double) 899 maxX(double) maxY(double) maxZ(double) 900 numPhysicalTags(size_t) physicalTag(int) ... 901 numBoundingPoints(size_t) pointTag(int) ... 902 ... 903 surfaceTag(int) minX(double) minY(double) minZ(double) 904 maxX(double) maxY(double) maxZ(double) 905 numPhysicalTags(size_t) physicalTag(int) ... 906 numBoundingCurves(size_t) curveTag(int) ... 907 ... 908 volumeTag(int) minX(double) minY(double) minZ(double) 909 maxX(double) maxY(double) maxZ(double) 910 numPhysicalTags(size_t) physicalTag(int) ... 911 numBoundngSurfaces(size_t) surfaceTag(int) ... 912 ... 913 $EndEntities 914 */ 915 static PetscErrorCode GmshReadEntities_v41(GmshFile *gmsh, GmshMesh *mesh) 916 { 917 PetscInt count[4], index, numTags, i; 918 int dim, eid, *tags = NULL; 919 GmshEntity *entity = NULL; 920 PetscErrorCode ierr; 921 922 PetscFunctionBegin; 923 ierr = GmshReadSize(gmsh, count, 4);CHKERRQ(ierr); 924 ierr = GmshEntitiesCreate(count, &mesh->entities);CHKERRQ(ierr); 925 for (dim = 0; dim < 4; ++dim) { 926 for (index = 0; index < count[dim]; ++index) { 927 ierr = GmshReadInt(gmsh, &eid, 1);CHKERRQ(ierr); 928 ierr = GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity);CHKERRQ(ierr); 929 ierr = GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6);CHKERRQ(ierr); 930 ierr = GmshReadSize(gmsh, &numTags, 1);CHKERRQ(ierr); 931 ierr = GmshBufferGet(gmsh, numTags, sizeof(int), &tags);CHKERRQ(ierr); 932 ierr = GmshReadInt(gmsh, tags, numTags);CHKERRQ(ierr); 933 entity->numTags = PetscMin(numTags, 4); 934 for (i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i]; 935 if (dim == 0) continue; 936 ierr = GmshReadSize(gmsh, &numTags, 1);CHKERRQ(ierr); 937 ierr = GmshBufferGet(gmsh, numTags, sizeof(int), &tags);CHKERRQ(ierr); 938 ierr = GmshReadInt(gmsh, tags, numTags);CHKERRQ(ierr); 939 } 940 } 941 PetscFunctionReturn(0); 942 } 943 944 /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format 945 $Nodes 946 numEntityBlocks(size_t) numNodes(size_t) 947 minNodeTag(size_t) maxNodeTag(size_t) 948 entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesBlock(size_t) 949 nodeTag(size_t) 950 ... 951 x(double) y(double) z(double) 952 < u(double; if parametric and entityDim = 1 or entityDim = 2) > 953 < v(double; if parametric and entityDim = 2) > 954 ... 955 ... 956 $EndNodes 957 */ 958 static PetscErrorCode GmshReadNodes_v41(GmshFile *gmsh, GmshMesh *mesh) 959 { 960 int info[3]; 961 PetscInt sizes[4], numEntityBlocks, numNodes, numNodesBlock = 0, block, node; 962 GmshNodes *nodes; 963 PetscErrorCode ierr; 964 965 PetscFunctionBegin; 966 ierr = GmshReadSize(gmsh, sizes, 4);CHKERRQ(ierr); 967 numEntityBlocks = sizes[0]; numNodes = sizes[1]; 968 ierr = GmshNodesCreate(numNodes, &nodes);CHKERRQ(ierr); 969 mesh->numNodes = numNodes; 970 mesh->nodelist = nodes; 971 for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) { 972 ierr = GmshReadInt(gmsh, info, 3);CHKERRQ(ierr); 973 PetscCheckFalse(info[2] != 0,PETSC_COMM_SELF, PETSC_ERR_SUP, "Parametric coordinates not supported"); 974 ierr = GmshReadSize(gmsh, &numNodesBlock, 1);CHKERRQ(ierr); 975 ierr = GmshReadSize(gmsh, nodes->id+node, numNodesBlock);CHKERRQ(ierr); 976 ierr = GmshReadDouble(gmsh, nodes->xyz+node*3, numNodesBlock*3);CHKERRQ(ierr); 977 } 978 gmsh->nodeStart = sizes[2]; 979 gmsh->nodeEnd = sizes[3]+1; 980 PetscFunctionReturn(0); 981 } 982 983 /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format 984 $Elements 985 numEntityBlocks(size_t) numElements(size_t) 986 minElementTag(size_t) maxElementTag(size_t) 987 entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t) 988 elementTag(size_t) nodeTag(size_t) ... 989 ... 990 ... 991 $EndElements 992 */ 993 static PetscErrorCode GmshReadElements_v41(GmshFile *gmsh, GmshMesh *mesh) 994 { 995 int info[3], eid, dim, cellType; 996 PetscInt sizes[4], *ibuf = NULL, numEntityBlocks, numElements, numBlockElements, numVerts, numNodes, numTags, block, elem, c, p; 997 GmshEntity *entity = NULL; 998 GmshElement *elements; 999 PetscInt *nodeMap = gmsh->nodeMap; 1000 PetscErrorCode ierr; 1001 1002 PetscFunctionBegin; 1003 ierr = GmshReadSize(gmsh, sizes, 4);CHKERRQ(ierr); 1004 numEntityBlocks = sizes[0]; numElements = sizes[1]; 1005 ierr = GmshElementsCreate(numElements, &elements);CHKERRQ(ierr); 1006 mesh->numElems = numElements; 1007 mesh->elements = elements; 1008 for (c = 0, block = 0; block < numEntityBlocks; ++block) { 1009 ierr = GmshReadInt(gmsh, info, 3);CHKERRQ(ierr); 1010 dim = info[0]; eid = info[1]; cellType = info[2]; 1011 ierr = GmshEntitiesGet(mesh->entities, dim, eid, &entity);CHKERRQ(ierr); 1012 ierr = GmshCellTypeCheck(cellType);CHKERRQ(ierr); 1013 numVerts = GmshCellMap[cellType].numVerts; 1014 numNodes = GmshCellMap[cellType].numNodes; 1015 numTags = entity->numTags; 1016 ierr = GmshReadSize(gmsh, &numBlockElements, 1);CHKERRQ(ierr); 1017 ierr = GmshBufferGet(gmsh, (1+numNodes)*numBlockElements, sizeof(PetscInt), &ibuf);CHKERRQ(ierr); 1018 ierr = GmshReadSize(gmsh, ibuf, (1+numNodes)*numBlockElements);CHKERRQ(ierr); 1019 for (elem = 0; elem < numBlockElements; ++elem, ++c) { 1020 GmshElement *element = elements + c; 1021 const PetscInt *id = ibuf + elem*(1+numNodes), *nodes = id + 1; 1022 element->id = id[0]; 1023 element->dim = dim; 1024 element->cellType = cellType; 1025 element->numVerts = numVerts; 1026 element->numNodes = numNodes; 1027 element->numTags = numTags; 1028 ierr = PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);CHKERRQ(ierr); 1029 for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]]; 1030 for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p]; 1031 } 1032 } 1033 PetscFunctionReturn(0); 1034 } 1035 1036 /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format 1037 $Periodic 1038 numPeriodicLinks(size_t) 1039 entityDim(int) entityTag(int) entityTagPrimary(int) 1040 numAffine(size_t) value(double) ... 1041 numCorrespondingNodes(size_t) 1042 nodeTag(size_t) nodeTagPrimary(size_t) 1043 ... 1044 ... 1045 $EndPeriodic 1046 */ 1047 static PetscErrorCode GmshReadPeriodic_v41(GmshFile *gmsh, PetscInt periodicMap[]) 1048 { 1049 int info[3]; 1050 double dbuf[16]; 1051 PetscInt numPeriodicLinks, numAffine, numCorrespondingNodes, *nodeTags = NULL, link, node; 1052 PetscInt *nodeMap = gmsh->nodeMap; 1053 PetscErrorCode ierr; 1054 1055 PetscFunctionBegin; 1056 ierr = GmshReadSize(gmsh, &numPeriodicLinks, 1);CHKERRQ(ierr); 1057 for (link = 0; link < numPeriodicLinks; ++link) { 1058 ierr = GmshReadInt(gmsh, info, 3);CHKERRQ(ierr); 1059 ierr = GmshReadSize(gmsh, &numAffine, 1);CHKERRQ(ierr); 1060 ierr = GmshReadDouble(gmsh, dbuf, numAffine);CHKERRQ(ierr); 1061 ierr = GmshReadSize(gmsh, &numCorrespondingNodes, 1);CHKERRQ(ierr); 1062 ierr = GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags);CHKERRQ(ierr); 1063 ierr = GmshReadSize(gmsh, nodeTags, numCorrespondingNodes*2);CHKERRQ(ierr); 1064 for (node = 0; node < numCorrespondingNodes; ++node) { 1065 PetscInt correspondingNode = nodeMap[nodeTags[node*2+0]]; 1066 PetscInt primaryNode = nodeMap[nodeTags[node*2+1]]; 1067 periodicMap[correspondingNode] = primaryNode; 1068 } 1069 } 1070 PetscFunctionReturn(0); 1071 } 1072 1073 /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format 1074 $MeshFormat // same as MSH version 2 1075 version(ASCII double; currently 4.1) 1076 file-type(ASCII int; 0 for ASCII mode, 1 for binary mode) 1077 data-size(ASCII int; sizeof(size_t)) 1078 < int with value one; only in binary mode, to detect endianness > 1079 $EndMeshFormat 1080 */ 1081 static PetscErrorCode GmshReadMeshFormat(GmshFile *gmsh) 1082 { 1083 char line[PETSC_MAX_PATH_LEN]; 1084 int snum, fileType, fileFormat, dataSize, checkEndian; 1085 float version; 1086 PetscErrorCode ierr; 1087 1088 PetscFunctionBegin; 1089 ierr = GmshReadString(gmsh, line, 3);CHKERRQ(ierr); 1090 snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize); 1091 PetscCheckFalse(snum != 3,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line); 1092 PetscCheckFalse(version < 2.2,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version); 1093 PetscCheckFalse((int)version == 3,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version); 1094 PetscCheckFalse(version > 4.1,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version); 1095 PetscCheckFalse(gmsh->binary && !fileType,PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is binary but Gmsh file is ASCII"); 1096 PetscCheckFalse(!gmsh->binary && fileType,PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is ASCII but Gmsh file is binary"); 1097 fileFormat = (int)roundf(version*10); 1098 PetscCheckFalse(fileFormat <= 40 && dataSize != sizeof(double),PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Data size %d is not valid for a Gmsh file", dataSize); 1099 PetscCheckFalse(fileFormat >= 41 && dataSize != sizeof(int) && dataSize != sizeof(PetscInt64),PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Data size %d is not valid for a Gmsh file", dataSize); 1100 gmsh->fileFormat = fileFormat; 1101 gmsh->dataSize = dataSize; 1102 gmsh->byteSwap = PETSC_FALSE; 1103 if (gmsh->binary) { 1104 ierr = GmshReadInt(gmsh, &checkEndian, 1);CHKERRQ(ierr); 1105 if (checkEndian != 1) { 1106 ierr = PetscByteSwap(&checkEndian, PETSC_ENUM, 1);CHKERRQ(ierr); 1107 PetscCheckFalse(checkEndian != 1,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to detect endianness in Gmsh file header: %s", line); 1108 gmsh->byteSwap = PETSC_TRUE; 1109 } 1110 } 1111 PetscFunctionReturn(0); 1112 } 1113 1114 /* 1115 PhysicalNames 1116 numPhysicalNames(ASCII int) 1117 dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max) 1118 ... 1119 $EndPhysicalNames 1120 */ 1121 static PetscErrorCode GmshReadPhysicalNames(GmshFile *gmsh, GmshMesh *mesh) 1122 { 1123 char line[PETSC_MAX_PATH_LEN], name[128+2], *p, *q; 1124 int snum, region, dim, tag; 1125 PetscErrorCode ierr; 1126 1127 PetscFunctionBegin; 1128 ierr = GmshReadString(gmsh, line, 1);CHKERRQ(ierr); 1129 snum = sscanf(line, "%d", ®ion); 1130 mesh->numRegions = region; 1131 PetscCheckFalse(snum != 1,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 1132 ierr = PetscMalloc2(mesh->numRegions, &mesh->regionTags, mesh->numRegions, &mesh->regionNames);CHKERRQ(ierr); 1133 for (region = 0; region < mesh->numRegions; ++region) { 1134 ierr = GmshReadString(gmsh, line, 2);CHKERRQ(ierr); 1135 snum = sscanf(line, "%d %d", &dim, &tag); 1136 PetscCheckFalse(snum != 2,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 1137 ierr = GmshReadString(gmsh, line, -(PetscInt)sizeof(line));CHKERRQ(ierr); 1138 ierr = PetscStrchr(line, '"', &p);CHKERRQ(ierr); 1139 PetscCheckFalse(!p,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 1140 ierr = PetscStrrchr(line, '"', &q);CHKERRQ(ierr); 1141 PetscCheckFalse(q == p,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 1142 ierr = PetscStrncpy(name, p+1, (size_t)(q-p-1));CHKERRQ(ierr); 1143 mesh->regionTags[region] = tag; 1144 ierr = PetscStrallocpy(name, &mesh->regionNames[region]);CHKERRQ(ierr); 1145 } 1146 PetscFunctionReturn(0); 1147 } 1148 1149 static PetscErrorCode GmshReadEntities(GmshFile *gmsh, GmshMesh *mesh) 1150 { 1151 PetscErrorCode ierr; 1152 1153 PetscFunctionBegin; 1154 switch (gmsh->fileFormat) { 1155 case 41: ierr = GmshReadEntities_v41(gmsh, mesh);CHKERRQ(ierr); break; 1156 default: ierr = GmshReadEntities_v40(gmsh, mesh);CHKERRQ(ierr); break; 1157 } 1158 PetscFunctionReturn(0); 1159 } 1160 1161 static PetscErrorCode GmshReadNodes(GmshFile *gmsh, GmshMesh *mesh) 1162 { 1163 PetscErrorCode ierr; 1164 1165 PetscFunctionBegin; 1166 switch (gmsh->fileFormat) { 1167 case 41: ierr = GmshReadNodes_v41(gmsh, mesh);CHKERRQ(ierr); break; 1168 case 40: ierr = GmshReadNodes_v40(gmsh, mesh);CHKERRQ(ierr); break; 1169 default: ierr = GmshReadNodes_v22(gmsh, mesh);CHKERRQ(ierr); break; 1170 } 1171 1172 { /* Gmsh v2.2/v4.0 does not provide min/max node tags */ 1173 if (mesh->numNodes > 0 && gmsh->nodeEnd >= gmsh->nodeStart) { 1174 PetscInt tagMin = PETSC_MAX_INT, tagMax = PETSC_MIN_INT, n; 1175 GmshNodes *nodes = mesh->nodelist; 1176 for (n = 0; n < mesh->numNodes; ++n) { 1177 const PetscInt tag = nodes->id[n]; 1178 tagMin = PetscMin(tag, tagMin); 1179 tagMax = PetscMax(tag, tagMax); 1180 } 1181 gmsh->nodeStart = tagMin; 1182 gmsh->nodeEnd = tagMax+1; 1183 } 1184 } 1185 1186 { /* Support for sparse node tags */ 1187 PetscInt n, t; 1188 GmshNodes *nodes = mesh->nodelist; 1189 ierr = PetscMalloc1(gmsh->nodeEnd - gmsh->nodeStart, &gmsh->nbuf);CHKERRQ(ierr); 1190 for (t = 0; t < gmsh->nodeEnd - gmsh->nodeStart; ++t) gmsh->nbuf[t] = PETSC_MIN_INT; 1191 gmsh->nodeMap = gmsh->nbuf - gmsh->nodeStart; 1192 for (n = 0; n < mesh->numNodes; ++n) { 1193 const PetscInt tag = nodes->id[n]; 1194 PetscCheckFalse(gmsh->nodeMap[tag] >= 0,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Repeated node tag %D", tag); 1195 gmsh->nodeMap[tag] = n; 1196 } 1197 } 1198 PetscFunctionReturn(0); 1199 } 1200 1201 static PetscErrorCode GmshReadElements(GmshFile *gmsh, GmshMesh *mesh) 1202 { 1203 PetscErrorCode ierr; 1204 1205 PetscFunctionBegin; 1206 switch (gmsh->fileFormat) { 1207 case 41: ierr = GmshReadElements_v41(gmsh, mesh);CHKERRQ(ierr); break; 1208 case 40: ierr = GmshReadElements_v40(gmsh, mesh);CHKERRQ(ierr); break; 1209 default: ierr = GmshReadElements_v22(gmsh, mesh);CHKERRQ(ierr); break; 1210 } 1211 1212 { /* Reorder elements by codimension and polytope type */ 1213 PetscInt ne = mesh->numElems; 1214 GmshElement *elements = mesh->elements; 1215 PetscInt keymap[GMSH_NUM_POLYTOPES], nk = 0; 1216 PetscInt offset[GMSH_NUM_POLYTOPES+1], e, k; 1217 1218 for (k = 0; k < GMSH_NUM_POLYTOPES; ++k) keymap[k] = PETSC_MIN_INT; 1219 ierr = PetscMemzero(offset,sizeof(offset));CHKERRQ(ierr); 1220 1221 keymap[GMSH_TET] = nk++; 1222 keymap[GMSH_HEX] = nk++; 1223 keymap[GMSH_PRI] = nk++; 1224 keymap[GMSH_PYR] = nk++; 1225 keymap[GMSH_TRI] = nk++; 1226 keymap[GMSH_QUA] = nk++; 1227 keymap[GMSH_SEG] = nk++; 1228 keymap[GMSH_VTX] = nk++; 1229 1230 ierr = GmshElementsCreate(mesh->numElems, &mesh->elements);CHKERRQ(ierr); 1231 #define key(eid) keymap[GmshCellMap[elements[(eid)].cellType].polytope] 1232 for (e = 0; e < ne; ++e) offset[1+key(e)]++; 1233 for (k = 1; k < nk; ++k) offset[k] += offset[k-1]; 1234 for (e = 0; e < ne; ++e) mesh->elements[offset[key(e)]++] = elements[e]; 1235 #undef key 1236 ierr = GmshElementsDestroy(&elements);CHKERRQ(ierr); 1237 } 1238 1239 { /* Mesh dimension and order */ 1240 GmshElement *elem = mesh->numElems ? mesh->elements : NULL; 1241 mesh->dim = elem ? GmshCellMap[elem->cellType].dim : 0; 1242 mesh->order = elem ? GmshCellMap[elem->cellType].order : 0; 1243 } 1244 1245 { 1246 PetscBT vtx; 1247 PetscInt dim = mesh->dim, e, n, v; 1248 1249 ierr = PetscBTCreate(mesh->numNodes, &vtx);CHKERRQ(ierr); 1250 1251 /* Compute number of cells and set of vertices */ 1252 mesh->numCells = 0; 1253 for (e = 0; e < mesh->numElems; ++e) { 1254 GmshElement *elem = mesh->elements + e; 1255 if (elem->dim == dim && dim > 0) mesh->numCells++; 1256 for (v = 0; v < elem->numVerts; v++) { 1257 ierr = PetscBTSet(vtx, elem->nodes[v]);CHKERRQ(ierr); 1258 } 1259 } 1260 1261 /* Compute numbering for vertices */ 1262 mesh->numVerts = 0; 1263 ierr = PetscMalloc1(mesh->numNodes, &mesh->vertexMap);CHKERRQ(ierr); 1264 for (n = 0; n < mesh->numNodes; ++n) 1265 mesh->vertexMap[n] = PetscBTLookup(vtx, n) ? mesh->numVerts++ : PETSC_MIN_INT; 1266 1267 ierr = PetscBTDestroy(&vtx);CHKERRQ(ierr); 1268 } 1269 PetscFunctionReturn(0); 1270 } 1271 1272 static PetscErrorCode GmshReadPeriodic(GmshFile *gmsh, GmshMesh *mesh) 1273 { 1274 PetscInt n; 1275 PetscErrorCode ierr; 1276 1277 PetscFunctionBegin; 1278 ierr = PetscMalloc1(mesh->numNodes, &mesh->periodMap);CHKERRQ(ierr); 1279 for (n = 0; n < mesh->numNodes; ++n) mesh->periodMap[n] = n; 1280 switch (gmsh->fileFormat) { 1281 case 41: ierr = GmshReadPeriodic_v41(gmsh, mesh->periodMap);CHKERRQ(ierr); break; 1282 default: ierr = GmshReadPeriodic_v40(gmsh, mesh->periodMap);CHKERRQ(ierr); break; 1283 } 1284 1285 /* Find canonical primary nodes */ 1286 for (n = 0; n < mesh->numNodes; ++n) 1287 while (mesh->periodMap[n] != mesh->periodMap[mesh->periodMap[n]]) 1288 mesh->periodMap[n] = mesh->periodMap[mesh->periodMap[n]]; 1289 1290 /* Renumber vertices (filter out correspondings) */ 1291 mesh->numVerts = 0; 1292 for (n = 0; n < mesh->numNodes; ++n) 1293 if (mesh->vertexMap[n] >= 0) /* is vertex */ 1294 if (mesh->periodMap[n] == n) /* is primary */ 1295 mesh->vertexMap[n] = mesh->numVerts++; 1296 for (n = 0; n < mesh->numNodes; ++n) 1297 if (mesh->vertexMap[n] >= 0) /* is vertex */ 1298 if (mesh->periodMap[n] != n) /* is corresponding */ 1299 mesh->vertexMap[n] = mesh->vertexMap[mesh->periodMap[n]]; 1300 PetscFunctionReturn(0); 1301 } 1302 1303 #define DM_POLYTOPE_VERTEX DM_POLYTOPE_POINT 1304 #define DM_POLYTOPE_PYRAMID DM_POLYTOPE_UNKNOWN 1305 static const DMPolytopeType DMPolytopeMap[] = { 1306 /* GMSH_VTX */ DM_POLYTOPE_VERTEX, 1307 /* GMSH_SEG */ DM_POLYTOPE_SEGMENT, 1308 /* GMSH_TRI */ DM_POLYTOPE_TRIANGLE, 1309 /* GMSH_QUA */ DM_POLYTOPE_QUADRILATERAL, 1310 /* GMSH_TET */ DM_POLYTOPE_TETRAHEDRON, 1311 /* GMSH_HEX */ DM_POLYTOPE_HEXAHEDRON, 1312 /* GMSH_PRI */ DM_POLYTOPE_TRI_PRISM, 1313 /* GMSH_PYR */ DM_POLYTOPE_PYRAMID, 1314 DM_POLYTOPE_UNKNOWN 1315 }; 1316 1317 static inline DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt cellType) 1318 { 1319 return DMPolytopeMap[GmshCellMap[cellType].polytope]; 1320 } 1321 1322 static PetscErrorCode GmshCreateFE(MPI_Comm comm, const char prefix[], PetscBool isSimplex, PetscBool continuity, PetscDTNodeType nodeType, PetscInt dim, PetscInt Nc, PetscInt k, PetscFE *fem) 1323 { 1324 DM K; 1325 PetscSpace P; 1326 PetscDualSpace Q; 1327 PetscQuadrature q, fq; 1328 PetscBool isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1329 PetscBool endpoint = PETSC_TRUE; 1330 char name[32]; 1331 PetscErrorCode ierr; 1332 1333 PetscFunctionBegin; 1334 /* Create space */ 1335 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 1336 ierr = PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 1337 ierr = PetscSpacePolynomialSetTensor(P, isTensor);CHKERRQ(ierr); 1338 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 1339 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 1340 ierr = PetscSpaceSetDegree(P, k, PETSC_DETERMINE);CHKERRQ(ierr); 1341 if (prefix) { 1342 ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix);CHKERRQ(ierr); 1343 ierr = PetscSpaceSetFromOptions(P);CHKERRQ(ierr); 1344 ierr = PetscObjectSetOptionsPrefix((PetscObject) P, NULL);CHKERRQ(ierr); 1345 ierr = PetscSpaceGetDegree(P, &k, NULL);CHKERRQ(ierr); 1346 } 1347 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 1348 /* Create dual space */ 1349 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 1350 ierr = PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1351 ierr = PetscDualSpaceLagrangeSetTensor(Q, isTensor);CHKERRQ(ierr); 1352 ierr = PetscDualSpaceLagrangeSetContinuity(Q, continuity);CHKERRQ(ierr); 1353 ierr = PetscDualSpaceLagrangeSetNodeType(Q, nodeType, endpoint, 0);CHKERRQ(ierr); 1354 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 1355 ierr = PetscDualSpaceSetOrder(Q, k);CHKERRQ(ierr); 1356 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 1357 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 1358 ierr = DMDestroy(&K);CHKERRQ(ierr); 1359 if (prefix) { 1360 ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);CHKERRQ(ierr); 1361 ierr = PetscDualSpaceSetFromOptions(Q);CHKERRQ(ierr); 1362 ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, NULL);CHKERRQ(ierr); 1363 } 1364 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 1365 /* Create quadrature */ 1366 if (isSimplex) { 1367 ierr = PetscDTStroudConicalQuadrature(dim, 1, k+1, -1, +1, &q);CHKERRQ(ierr); 1368 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, k+1, -1, +1, &fq);CHKERRQ(ierr); 1369 } else { 1370 ierr = PetscDTGaussTensorQuadrature(dim, 1, k+1, -1, +1, &q);CHKERRQ(ierr); 1371 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, k+1, -1, +1, &fq);CHKERRQ(ierr); 1372 } 1373 /* Create finite element */ 1374 ierr = PetscFECreate(comm, fem);CHKERRQ(ierr); 1375 ierr = PetscFESetType(*fem, PETSCFEBASIC);CHKERRQ(ierr); 1376 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 1377 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 1378 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 1379 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 1380 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 1381 if (prefix) { 1382 ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);CHKERRQ(ierr); 1383 ierr = PetscFESetFromOptions(*fem);CHKERRQ(ierr); 1384 ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, NULL);CHKERRQ(ierr); 1385 } 1386 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 1387 /* Cleanup */ 1388 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 1389 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 1390 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 1391 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 1392 /* Set finite element name */ 1393 ierr = PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);CHKERRQ(ierr); 1394 ierr = PetscFESetName(*fem, name);CHKERRQ(ierr); 1395 PetscFunctionReturn(0); 1396 } 1397 1398 /*@C 1399 DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file 1400 1401 + comm - The MPI communicator 1402 . filename - Name of the Gmsh file 1403 - interpolate - Create faces and edges in the mesh 1404 1405 Output Parameter: 1406 . dm - The DM object representing the mesh 1407 1408 Level: beginner 1409 1410 .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 1411 @*/ 1412 PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 1413 { 1414 PetscViewer viewer; 1415 PetscMPIInt rank; 1416 int fileType; 1417 PetscViewerType vtype; 1418 PetscErrorCode ierr; 1419 1420 PetscFunctionBegin; 1421 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 1422 1423 /* Determine Gmsh file type (ASCII or binary) from file header */ 1424 if (rank == 0) { 1425 GmshFile gmsh[1]; 1426 char line[PETSC_MAX_PATH_LEN]; 1427 int snum; 1428 float version; 1429 1430 ierr = PetscArrayzero(gmsh,1);CHKERRQ(ierr); 1431 ierr = PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer);CHKERRQ(ierr); 1432 ierr = PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 1433 ierr = PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ);CHKERRQ(ierr); 1434 ierr = PetscViewerFileSetName(gmsh->viewer, filename);CHKERRQ(ierr); 1435 /* Read only the first two lines of the Gmsh file */ 1436 ierr = GmshReadSection(gmsh, line);CHKERRQ(ierr); 1437 ierr = GmshExpect(gmsh, "$MeshFormat", line);CHKERRQ(ierr); 1438 ierr = GmshReadString(gmsh, line, 2);CHKERRQ(ierr); 1439 snum = sscanf(line, "%f %d", &version, &fileType); 1440 PetscCheckFalse(snum != 2,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line); 1441 PetscCheckFalse(version < 2.2,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version); 1442 PetscCheckFalse((int)version == 3,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version); 1443 PetscCheckFalse(version > 4.1,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version); 1444 ierr = PetscViewerDestroy(&gmsh->viewer);CHKERRQ(ierr); 1445 } 1446 ierr = MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);CHKERRMPI(ierr); 1447 vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY; 1448 1449 /* Create appropriate viewer and build plex */ 1450 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr); 1451 ierr = PetscViewerSetType(viewer, vtype);CHKERRQ(ierr); 1452 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr); 1453 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr); 1454 ierr = DMPlexCreateGmsh(comm, viewer, interpolate, dm);CHKERRQ(ierr); 1455 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1456 PetscFunctionReturn(0); 1457 } 1458 1459 /*@ 1460 DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer 1461 1462 Collective 1463 1464 Input Parameters: 1465 + comm - The MPI communicator 1466 . viewer - The Viewer associated with a Gmsh file 1467 - interpolate - Create faces and edges in the mesh 1468 1469 Output Parameter: 1470 . dm - The DM object representing the mesh 1471 1472 Note: http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format 1473 1474 Level: beginner 1475 1476 .seealso: DMPLEX, DMCreate() 1477 @*/ 1478 PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm) 1479 { 1480 GmshMesh *mesh = NULL; 1481 PetscViewer parentviewer = NULL; 1482 PetscBT periodicVerts = NULL; 1483 PetscBT periodicCells = NULL; 1484 DM cdm; 1485 PetscSection coordSection; 1486 Vec coordinates; 1487 DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL, marker = NULL, *regionSets; 1488 PetscInt dim = 0, coordDim = -1, order = 0; 1489 PetscInt numNodes = 0, numElems = 0, numVerts = 0, numCells = 0; 1490 PetscInt cell, cone[8], e, n, v, d; 1491 PetscBool binary, usemarker = PETSC_FALSE, useregions = PETSC_FALSE; 1492 PetscBool hybrid = interpolate, periodic = PETSC_TRUE; 1493 PetscBool highOrder = PETSC_TRUE, highOrderSet, project = PETSC_FALSE; 1494 PetscBool isSimplex = PETSC_FALSE, isHybrid = PETSC_FALSE, hasTetra = PETSC_FALSE; 1495 PetscMPIInt rank; 1496 PetscErrorCode ierr; 1497 1498 PetscFunctionBegin; 1499 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 1500 ierr = PetscObjectOptionsBegin((PetscObject)viewer);CHKERRQ(ierr); 1501 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Gmsh options");CHKERRQ(ierr); 1502 ierr = PetscOptionsBool("-dm_plex_gmsh_hybrid", "Generate hybrid cell bounds", "DMPlexCreateGmsh", hybrid, &hybrid, NULL);CHKERRQ(ierr); 1503 ierr = PetscOptionsBool("-dm_plex_gmsh_periodic","Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL);CHKERRQ(ierr); 1504 ierr = PetscOptionsBool("-dm_plex_gmsh_highorder","Generate high-order coordinates", "DMPlexCreateGmsh", highOrder, &highOrder, &highOrderSet);CHKERRQ(ierr); 1505 ierr = PetscOptionsBool("-dm_plex_gmsh_project", "Project high-order coordinates to a different space", "DMPlexCreateGmsh", project, &project, NULL);CHKERRQ(ierr); 1506 ierr = PetscOptionsBool("-dm_plex_gmsh_use_marker", "Generate marker label", "DMPlexCreateGmsh", usemarker, &usemarker, NULL);CHKERRQ(ierr); 1507 ierr = PetscOptionsBool("-dm_plex_gmsh_use_regions", "Generate labels with region names", "DMPlexCreateGmsh", useregions, &useregions, NULL);CHKERRQ(ierr); 1508 ierr = PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE);CHKERRQ(ierr); 1509 ierr = PetscOptionsTail();CHKERRQ(ierr); 1510 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1511 1512 ierr = GmshCellInfoSetUp();CHKERRQ(ierr); 1513 1514 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 1515 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 1516 ierr = PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,NULL,NULL,NULL);CHKERRQ(ierr); 1517 1518 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);CHKERRQ(ierr); 1519 1520 /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */ 1521 if (binary) { 1522 parentviewer = viewer; 1523 ierr = PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);CHKERRQ(ierr); 1524 } 1525 1526 if (rank == 0) { 1527 GmshFile gmsh[1]; 1528 char line[PETSC_MAX_PATH_LEN]; 1529 PetscBool match; 1530 1531 ierr = PetscArrayzero(gmsh,1);CHKERRQ(ierr); 1532 gmsh->viewer = viewer; 1533 gmsh->binary = binary; 1534 1535 ierr = GmshMeshCreate(&mesh);CHKERRQ(ierr); 1536 1537 /* Read mesh format */ 1538 ierr = GmshReadSection(gmsh, line);CHKERRQ(ierr); 1539 ierr = GmshExpect(gmsh, "$MeshFormat", line);CHKERRQ(ierr); 1540 ierr = GmshReadMeshFormat(gmsh);CHKERRQ(ierr); 1541 ierr = GmshReadEndSection(gmsh, "$EndMeshFormat", line);CHKERRQ(ierr); 1542 1543 /* OPTIONAL Read physical names */ 1544 ierr = GmshReadSection(gmsh, line);CHKERRQ(ierr); 1545 ierr = GmshMatch(gmsh, "$PhysicalNames", line, &match);CHKERRQ(ierr); 1546 if (match) { 1547 ierr = GmshExpect(gmsh, "$PhysicalNames", line);CHKERRQ(ierr); 1548 ierr = GmshReadPhysicalNames(gmsh, mesh);CHKERRQ(ierr); 1549 ierr = GmshReadEndSection(gmsh, "$EndPhysicalNames", line);CHKERRQ(ierr); 1550 /* Initial read for entity section */ 1551 ierr = GmshReadSection(gmsh, line);CHKERRQ(ierr); 1552 } 1553 1554 /* Read entities */ 1555 if (gmsh->fileFormat >= 40) { 1556 ierr = GmshExpect(gmsh, "$Entities", line);CHKERRQ(ierr); 1557 ierr = GmshReadEntities(gmsh, mesh);CHKERRQ(ierr); 1558 ierr = GmshReadEndSection(gmsh, "$EndEntities", line);CHKERRQ(ierr); 1559 /* Initial read for nodes section */ 1560 ierr = GmshReadSection(gmsh, line);CHKERRQ(ierr); 1561 } 1562 1563 /* Read nodes */ 1564 ierr = GmshExpect(gmsh, "$Nodes", line);CHKERRQ(ierr); 1565 ierr = GmshReadNodes(gmsh, mesh);CHKERRQ(ierr); 1566 ierr = GmshReadEndSection(gmsh, "$EndNodes", line);CHKERRQ(ierr); 1567 1568 /* Read elements */ 1569 ierr = GmshReadSection(gmsh, line);CHKERRQ(ierr); 1570 ierr = GmshExpect(gmsh, "$Elements", line);CHKERRQ(ierr); 1571 ierr = GmshReadElements(gmsh, mesh);CHKERRQ(ierr); 1572 ierr = GmshReadEndSection(gmsh, "$EndElements", line);CHKERRQ(ierr); 1573 1574 /* Read periodic section (OPTIONAL) */ 1575 if (periodic) { 1576 ierr = GmshReadSection(gmsh, line);CHKERRQ(ierr); 1577 ierr = GmshMatch(gmsh, "$Periodic", line, &periodic);CHKERRQ(ierr); 1578 } 1579 if (periodic) { 1580 ierr = GmshExpect(gmsh, "$Periodic", line);CHKERRQ(ierr); 1581 ierr = GmshReadPeriodic(gmsh, mesh);CHKERRQ(ierr); 1582 ierr = GmshReadEndSection(gmsh, "$EndPeriodic", line);CHKERRQ(ierr); 1583 } 1584 1585 ierr = PetscFree(gmsh->wbuf);CHKERRQ(ierr); 1586 ierr = PetscFree(gmsh->sbuf);CHKERRQ(ierr); 1587 ierr = PetscFree(gmsh->nbuf);CHKERRQ(ierr); 1588 1589 dim = mesh->dim; 1590 order = mesh->order; 1591 numNodes = mesh->numNodes; 1592 numElems = mesh->numElems; 1593 numVerts = mesh->numVerts; 1594 numCells = mesh->numCells; 1595 1596 { 1597 GmshElement *elemA = mesh->numCells > 0 ? mesh->elements : NULL; 1598 GmshElement *elemB = elemA ? elemA + mesh->numCells - 1 : NULL; 1599 int ptA = elemA ? GmshCellMap[elemA->cellType].polytope : -1; 1600 int ptB = elemB ? GmshCellMap[elemB->cellType].polytope : -1; 1601 isSimplex = (ptA == GMSH_QUA || ptA == GMSH_HEX) ? PETSC_FALSE : PETSC_TRUE; 1602 isHybrid = (ptA == ptB) ? PETSC_FALSE : PETSC_TRUE; 1603 hasTetra = (ptA == GMSH_TET) ? PETSC_TRUE : PETSC_FALSE; 1604 } 1605 } 1606 1607 if (parentviewer) { 1608 ierr = PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);CHKERRQ(ierr); 1609 } 1610 1611 { 1612 int buf[6]; 1613 1614 buf[0] = (int)dim; 1615 buf[1] = (int)order; 1616 buf[2] = periodic; 1617 buf[3] = isSimplex; 1618 buf[4] = isHybrid; 1619 buf[5] = hasTetra; 1620 1621 ierr = MPI_Bcast(buf, 6, MPI_INT, 0, comm);CHKERRMPI(ierr); 1622 1623 dim = buf[0]; 1624 order = buf[1]; 1625 periodic = buf[2] ? PETSC_TRUE : PETSC_FALSE; 1626 isSimplex = buf[3] ? PETSC_TRUE : PETSC_FALSE; 1627 isHybrid = buf[4] ? PETSC_TRUE : PETSC_FALSE; 1628 hasTetra = buf[5] ? PETSC_TRUE : PETSC_FALSE; 1629 } 1630 1631 if (!highOrderSet) highOrder = (order > 1) ? PETSC_TRUE : PETSC_FALSE; 1632 PetscCheckFalse(highOrder && isHybrid,comm, PETSC_ERR_SUP, "No support for discretization on hybrid meshes yet"); 1633 1634 /* We do not want this label automatically computed, instead we fill it here */ 1635 ierr = DMCreateLabel(*dm, "celltype");CHKERRQ(ierr); 1636 1637 /* Allocate the cell-vertex mesh */ 1638 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr); 1639 for (cell = 0; cell < numCells; ++cell) { 1640 GmshElement *elem = mesh->elements + cell; 1641 DMPolytopeType ctype = DMPolytopeTypeFromGmsh(elem->cellType); 1642 if (hybrid && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR; 1643 ierr = DMPlexSetConeSize(*dm, cell, elem->numVerts);CHKERRQ(ierr); 1644 ierr = DMPlexSetCellType(*dm, cell, ctype);CHKERRQ(ierr); 1645 } 1646 for (v = numCells; v < numCells+numVerts; ++v) { 1647 ierr = DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr); 1648 } 1649 ierr = DMSetUp(*dm);CHKERRQ(ierr); 1650 1651 /* Add cell-vertex connections */ 1652 for (cell = 0; cell < numCells; ++cell) { 1653 GmshElement *elem = mesh->elements + cell; 1654 for (v = 0; v < elem->numVerts; ++v) { 1655 const PetscInt nn = elem->nodes[v]; 1656 const PetscInt vv = mesh->vertexMap[nn]; 1657 cone[v] = numCells + vv; 1658 } 1659 ierr = DMPlexReorderCell(*dm, cell, cone);CHKERRQ(ierr); 1660 ierr = DMPlexSetCone(*dm, cell, cone);CHKERRQ(ierr); 1661 } 1662 1663 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 1664 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 1665 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 1666 if (interpolate) { 1667 DM idm; 1668 1669 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 1670 ierr = DMDestroy(dm);CHKERRQ(ierr); 1671 *dm = idm; 1672 } 1673 1674 PetscCheckFalse(usemarker && !interpolate && dim > 1,comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation"); 1675 if (rank == 0 && usemarker) { 1676 PetscInt f, fStart, fEnd; 1677 1678 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr); 1679 ierr = DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1680 for (f = fStart; f < fEnd; ++f) { 1681 PetscInt suppSize; 1682 1683 ierr = DMPlexGetSupportSize(*dm, f, &suppSize);CHKERRQ(ierr); 1684 if (suppSize == 1) { 1685 PetscInt *cone = NULL, coneSize, p; 1686 1687 ierr = DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);CHKERRQ(ierr); 1688 for (p = 0; p < coneSize; p += 2) { 1689 ierr = DMSetLabelValue_Fast(*dm, &marker, "marker", cone[p], 1);CHKERRQ(ierr); 1690 } 1691 ierr = DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);CHKERRQ(ierr); 1692 } 1693 } 1694 } 1695 1696 if (rank == 0) { 1697 const PetscInt Nr = useregions ? mesh->numRegions : 0; 1698 PetscInt vStart, vEnd; 1699 1700 ierr = PetscCalloc1(Nr, ®ionSets);CHKERRQ(ierr); 1701 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 1702 for (cell = 0, e = 0; e < numElems; ++e) { 1703 GmshElement *elem = mesh->elements + e; 1704 1705 /* Create cell sets */ 1706 if (elem->dim == dim && dim > 0) { 1707 if (elem->numTags > 0) { 1708 const PetscInt tag = elem->tags[0]; 1709 PetscInt r; 1710 1711 ierr = DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", cell, tag);CHKERRQ(ierr); 1712 for (r = 0; r < Nr; ++r) { 1713 if (mesh->regionTags[r] == tag) {ierr = DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], cell, tag);CHKERRQ(ierr);} 1714 } 1715 } 1716 cell++; 1717 } 1718 1719 /* Create face sets */ 1720 if (interpolate && elem->dim == dim-1) { 1721 PetscInt joinSize; 1722 const PetscInt *join = NULL; 1723 const PetscInt tag = elem->tags[0]; 1724 PetscInt r; 1725 1726 /* Find the relevant facet with vertex joins */ 1727 for (v = 0; v < elem->numVerts; ++v) { 1728 const PetscInt nn = elem->nodes[v]; 1729 const PetscInt vv = mesh->vertexMap[nn]; 1730 cone[v] = vStart + vv; 1731 } 1732 ierr = DMPlexGetFullJoin(*dm, elem->numVerts, cone, &joinSize, &join);CHKERRQ(ierr); 1733 PetscCheckFalse(joinSize != 1,PETSC_COMM_SELF, PETSC_ERR_SUP, "Could not determine Plex facet for Gmsh element %D (Plex cell %D)", elem->id, e); 1734 ierr = DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], tag);CHKERRQ(ierr); 1735 for (r = 0; r < Nr; ++r) { 1736 if (mesh->regionTags[r] == tag) {ierr = DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], join[0], tag);CHKERRQ(ierr);} 1737 } 1738 ierr = DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join);CHKERRQ(ierr); 1739 } 1740 1741 /* Create vertex sets */ 1742 if (elem->dim == 0) { 1743 if (elem->numTags > 0) { 1744 const PetscInt nn = elem->nodes[0]; 1745 const PetscInt vv = mesh->vertexMap[nn]; 1746 const PetscInt tag = elem->tags[0]; 1747 PetscInt r; 1748 1749 ierr = DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag);CHKERRQ(ierr); 1750 for (r = 0; r < Nr; ++r) { 1751 if (mesh->regionTags[r] == tag) {ierr = DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag);CHKERRQ(ierr);} 1752 } 1753 } 1754 } 1755 } 1756 ierr = PetscFree(regionSets);CHKERRQ(ierr); 1757 } 1758 1759 { /* Create Cell/Face/Vertex Sets labels at all processes */ 1760 enum {n = 4}; 1761 PetscBool flag[n]; 1762 1763 flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE; 1764 flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE; 1765 flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE; 1766 flag[3] = marker ? PETSC_TRUE : PETSC_FALSE; 1767 ierr = MPI_Bcast(flag, n, MPIU_BOOL, 0, comm);CHKERRMPI(ierr); 1768 if (flag[0]) {ierr = DMCreateLabel(*dm, "Cell Sets");CHKERRQ(ierr);} 1769 if (flag[1]) {ierr = DMCreateLabel(*dm, "Face Sets");CHKERRQ(ierr);} 1770 if (flag[2]) {ierr = DMCreateLabel(*dm, "Vertex Sets");CHKERRQ(ierr);} 1771 if (flag[3]) {ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr);} 1772 } 1773 1774 if (periodic) { 1775 ierr = PetscBTCreate(numVerts, &periodicVerts);CHKERRQ(ierr); 1776 for (n = 0; n < numNodes; ++n) { 1777 if (mesh->vertexMap[n] >= 0) { 1778 if (PetscUnlikely(mesh->periodMap[n] != n)) { 1779 PetscInt m = mesh->periodMap[n]; 1780 ierr = PetscBTSet(periodicVerts, mesh->vertexMap[n]);CHKERRQ(ierr); 1781 ierr = PetscBTSet(periodicVerts, mesh->vertexMap[m]);CHKERRQ(ierr); 1782 } 1783 } 1784 } 1785 ierr = PetscBTCreate(numCells, &periodicCells);CHKERRQ(ierr); 1786 for (cell = 0; cell < numCells; ++cell) { 1787 GmshElement *elem = mesh->elements + cell; 1788 for (v = 0; v < elem->numVerts; ++v) { 1789 PetscInt nn = elem->nodes[v]; 1790 PetscInt vv = mesh->vertexMap[nn]; 1791 if (PetscUnlikely(PetscBTLookup(periodicVerts, vv))) { 1792 ierr = PetscBTSet(periodicCells, cell);CHKERRQ(ierr); break; 1793 } 1794 } 1795 } 1796 } 1797 1798 /* Setup coordinate DM */ 1799 if (coordDim < 0) coordDim = dim; 1800 ierr = DMSetCoordinateDim(*dm, coordDim);CHKERRQ(ierr); 1801 ierr = DMGetCoordinateDM(*dm, &cdm);CHKERRQ(ierr); 1802 if (highOrder) { 1803 PetscFE fe; 1804 PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE; 1805 PetscDTNodeType nodeType = PETSCDTNODES_EQUISPACED; 1806 1807 if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */ 1808 1809 ierr = GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);CHKERRQ(ierr); 1810 ierr = PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view");CHKERRQ(ierr); 1811 ierr = DMSetField(cdm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 1812 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 1813 ierr = DMCreateDS(cdm);CHKERRQ(ierr); 1814 } 1815 1816 /* Create coordinates */ 1817 if (highOrder) { 1818 1819 PetscInt maxDof = GmshNumNodes_HEX(order)*coordDim; 1820 double *coords = mesh ? mesh->nodelist->xyz : NULL; 1821 PetscSection section; 1822 PetscScalar *cellCoords; 1823 1824 ierr = DMSetLocalSection(cdm, NULL);CHKERRQ(ierr); 1825 ierr = DMGetLocalSection(cdm, &coordSection);CHKERRQ(ierr); 1826 ierr = PetscSectionClone(coordSection, §ion);CHKERRQ(ierr); 1827 ierr = DMPlexSetClosurePermutationTensor(cdm, 0, section);CHKERRQ(ierr); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */ 1828 1829 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 1830 ierr = PetscMalloc1(maxDof, &cellCoords);CHKERRQ(ierr); 1831 for (cell = 0; cell < numCells; ++cell) { 1832 GmshElement *elem = mesh->elements + cell; 1833 const int *lexorder = GmshCellMap[elem->cellType].lexorder(); 1834 int s = 0; 1835 for (n = 0; n < elem->numNodes; ++n) { 1836 while (lexorder[n+s] < 0) ++s; 1837 const PetscInt node = elem->nodes[lexorder[n+s]]; 1838 for (d = 0; d < coordDim; ++d) cellCoords[(n+s)*coordDim+d] = (PetscReal) coords[node*3+d]; 1839 } 1840 if (s) { 1841 /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */ 1842 PetscReal quaCenterWeights[9] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25}; 1843 /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */ 1844 PetscReal hexBottomWeights[27] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25, 1845 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1846 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; 1847 PetscReal hexFrontWeights[27] = {-0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1848 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1849 -0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; 1850 PetscReal hexLeftWeights[27] = {-0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 1851 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 1852 -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0}; 1853 PetscReal hexRightWeights[27] = { 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 1854 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1855 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25}; 1856 PetscReal hexBackWeights[27] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25, 1857 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 1858 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25}; 1859 PetscReal hexTopWeights[27] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1860 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1861 -0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25}; 1862 PetscReal hexCenterWeights[27] = {-0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25, 1863 0.25, 0.0, 0.25, 0.0, 0.0, 0.0, 0.25, 0.0, 0.25, 1864 -0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25}; 1865 PetscReal *sdWeights2[9] = {NULL, NULL, NULL, NULL, quaCenterWeights, NULL, NULL, NULL, NULL}; 1866 PetscReal *sdWeights3[27] = {NULL, NULL, NULL, NULL, hexBottomWeights, NULL, NULL, NULL, NULL, 1867 NULL, hexFrontWeights, NULL, hexLeftWeights, hexCenterWeights, hexRightWeights, NULL, hexBackWeights, NULL, 1868 NULL, NULL, NULL, NULL, hexTopWeights, NULL, NULL, NULL, NULL}; 1869 PetscReal **sdWeights[4] = {NULL, NULL, sdWeights2, sdWeights3}; 1870 1871 /* Missing entries in serendipity cell, only works for 8-node quad and 20-node hex */ 1872 for (n = 0; n < elem->numNodes+s; ++n) { 1873 if (lexorder[n] >= 0) continue; 1874 for (d = 0; d < coordDim; ++d) cellCoords[n*coordDim+d] = 0.0; 1875 for (int bn = 0; bn < elem->numNodes+s; ++bn) { 1876 if (lexorder[bn] < 0) continue; 1877 const PetscReal *weights = sdWeights[coordDim][n]; 1878 const PetscInt bnode = elem->nodes[lexorder[bn]]; 1879 for (d = 0; d < coordDim; ++d) cellCoords[n*coordDim+d] += weights[bn] * (PetscReal) coords[bnode*3+d]; 1880 } 1881 } 1882 } 1883 ierr = DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES);CHKERRQ(ierr); 1884 } 1885 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 1886 ierr = PetscFree(cellCoords);CHKERRQ(ierr); 1887 1888 } else { 1889 1890 PetscInt *nodeMap; 1891 double *coords = mesh ? mesh->nodelist->xyz : NULL; 1892 PetscScalar *pointCoords; 1893 1894 ierr = DMGetLocalSection(cdm, &coordSection);CHKERRQ(ierr); 1895 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 1896 ierr = PetscSectionSetFieldComponents(coordSection, 0, coordDim);CHKERRQ(ierr); 1897 if (periodic) { /* we need to localize coordinates on cells */ 1898 ierr = PetscSectionSetChart(coordSection, 0, numCells+numVerts);CHKERRQ(ierr); 1899 } else { 1900 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVerts);CHKERRQ(ierr); 1901 } 1902 if (periodic) { 1903 for (cell = 0; cell < numCells; ++cell) { 1904 if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) { 1905 GmshElement *elem = mesh->elements + cell; 1906 PetscInt dof = elem->numVerts * coordDim; 1907 ierr = PetscSectionSetDof(coordSection, cell, dof);CHKERRQ(ierr); 1908 ierr = PetscSectionSetFieldDof(coordSection, cell, 0, dof);CHKERRQ(ierr); 1909 } 1910 } 1911 } 1912 for (v = numCells; v < numCells+numVerts; ++v) { 1913 ierr = PetscSectionSetDof(coordSection, v, coordDim);CHKERRQ(ierr); 1914 ierr = PetscSectionSetFieldDof(coordSection, v, 0, coordDim);CHKERRQ(ierr); 1915 } 1916 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 1917 1918 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr); 1919 ierr = VecGetArray(coordinates, &pointCoords);CHKERRQ(ierr); 1920 if (periodic) { 1921 for (cell = 0; cell < numCells; ++cell) { 1922 if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) { 1923 GmshElement *elem = mesh->elements + cell; 1924 PetscInt off, node; 1925 for (v = 0; v < elem->numVerts; ++v) 1926 cone[v] = elem->nodes[v]; 1927 ierr = DMPlexReorderCell(cdm, cell, cone);CHKERRQ(ierr); 1928 ierr = PetscSectionGetOffset(coordSection, cell, &off);CHKERRQ(ierr); 1929 for (v = 0; v < elem->numVerts; ++v) 1930 for (node = cone[v], d = 0; d < coordDim; ++d) 1931 pointCoords[off++] = (PetscReal) coords[node*3+d]; 1932 } 1933 } 1934 } 1935 ierr = PetscMalloc1(numVerts, &nodeMap);CHKERRQ(ierr); 1936 for (n = 0; n < numNodes; n++) 1937 if (mesh->vertexMap[n] >= 0) 1938 nodeMap[mesh->vertexMap[n]] = n; 1939 for (v = 0; v < numVerts; ++v) { 1940 PetscInt off, node = nodeMap[v]; 1941 ierr = PetscSectionGetOffset(coordSection, numCells + v, &off);CHKERRQ(ierr); 1942 for (d = 0; d < coordDim; ++d) 1943 pointCoords[off+d] = (PetscReal) coords[node*3+d]; 1944 } 1945 ierr = PetscFree(nodeMap);CHKERRQ(ierr); 1946 ierr = VecRestoreArray(coordinates, &pointCoords);CHKERRQ(ierr); 1947 1948 } 1949 1950 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 1951 ierr = VecSetBlockSize(coordinates, coordDim);CHKERRQ(ierr); 1952 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 1953 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 1954 ierr = DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);CHKERRQ(ierr); 1955 1956 ierr = GmshMeshDestroy(&mesh);CHKERRQ(ierr); 1957 ierr = PetscBTDestroy(&periodicVerts);CHKERRQ(ierr); 1958 ierr = PetscBTDestroy(&periodicCells);CHKERRQ(ierr); 1959 1960 if (highOrder && project) { 1961 PetscFE fe; 1962 const char prefix[] = "dm_plex_gmsh_project_"; 1963 PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE; 1964 PetscDTNodeType nodeType = PETSCDTNODES_GAUSSJACOBI; 1965 1966 if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */ 1967 1968 ierr = GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);CHKERRQ(ierr); 1969 ierr = PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view");CHKERRQ(ierr); 1970 ierr = DMProjectCoordinates(*dm, fe);CHKERRQ(ierr); 1971 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 1972 } 1973 1974 ierr = PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,NULL,NULL,NULL);CHKERRQ(ierr); 1975 PetscFunctionReturn(0); 1976 } 1977