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