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