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_multiple_tags - Allow multiple tags for default labels 1559 - -dm_plex_gmsh_spacedim <d> - Embedding space dimension, if different from topological dimension 1560 1561 Level: beginner 1562 1563 Notes: 1564 The Gmsh file format is described in <http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format> 1565 1566 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`. 1567 1568 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()` 1569 @*/ 1570 PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm) 1571 { 1572 GmshMesh *mesh = NULL; 1573 PetscViewer parentviewer = NULL; 1574 PetscBT periodicVerts = NULL; 1575 PetscBT *periodicCells = NULL; 1576 DM cdm, cdmCell = NULL; 1577 PetscSection cs, csCell = NULL; 1578 Vec coordinates, coordinatesCell; 1579 DMLabel cellSets = NULL, faceSets = NULL, edgeSets = NULL, vertSets = NULL, marker = NULL, *regionSets; 1580 PetscInt dim = 0, coordDim = -1, order = 0, maxHeight = 0; 1581 PetscInt numNodes = 0, numElems = 0, numVerts = 0, numCells = 0, vStart, vEnd; 1582 PetscInt cell, cone[8], e, n, v, d; 1583 PetscBool usegeneric = PETSC_TRUE, useregions = PETSC_FALSE, markvertices = PETSC_FALSE, multipleTags = PETSC_FALSE; 1584 PetscBool flg, binary, hybrid = interpolate, periodic = PETSC_TRUE; 1585 PetscBool highOrder = PETSC_TRUE, highOrderSet, project = PETSC_FALSE; 1586 PetscBool isSimplex = PETSC_FALSE, isHybrid = PETSC_FALSE, hasTetra = PETSC_FALSE; 1587 PetscMPIInt rank; 1588 1589 PetscFunctionBegin; 1590 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1591 PetscObjectOptionsBegin((PetscObject)viewer); 1592 PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Gmsh options"); 1593 PetscCall(PetscOptionsBool("-dm_plex_gmsh_hybrid", "Force triangular prisms to use tensor order", "DMPlexCreateGmsh", hybrid, &hybrid, NULL)); 1594 PetscCall(PetscOptionsBool("-dm_plex_gmsh_periodic", "Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL)); 1595 PetscCall(PetscOptionsBool("-dm_plex_gmsh_highorder", "Generate high-order coordinates", "DMPlexCreateGmsh", highOrder, &highOrder, &highOrderSet)); 1596 PetscCall(PetscOptionsBool("-dm_plex_gmsh_project", "Project high-order coordinates to a different space", "DMPlexCreateGmsh", project, &project, NULL)); 1597 PetscCall(PetscOptionsBool("-dm_plex_gmsh_use_regions", "Generate labels with region names", "DMPlexCreateGmsh", useregions, &useregions, NULL)); 1598 PetscCall(PetscOptionsBool("-dm_plex_gmsh_use_generic", "Generate generic labels, i.e. Cell Sets, Face Sets, etc", "DMPlexCreateGmsh", usegeneric, &usegeneric, &flg)); 1599 if (!flg && useregions) usegeneric = PETSC_FALSE; 1600 PetscCall(PetscOptionsBool("-dm_plex_gmsh_mark_vertices", "Add vertices to generated labels", "DMPlexCreateGmsh", markvertices, &markvertices, NULL)); 1601 PetscCall(PetscOptionsBool("-dm_plex_gmsh_multiple_tags", "Allow multiple tags for default labels", "DMPlexCreateGmsh", multipleTags, &multipleTags, NULL)); 1602 PetscCall(PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE)); 1603 PetscCall(PetscOptionsBoundedInt("-dm_localize_height", "Localize edges and faces in addition to cells", "", maxHeight, &maxHeight, NULL, 0)); 1604 PetscOptionsHeadEnd(); 1605 PetscOptionsEnd(); 1606 1607 PetscCall(GmshCellInfoSetUp()); 1608 1609 PetscCall(DMCreate(comm, dm)); 1610 PetscCall(DMSetType(*dm, DMPLEX)); 1611 PetscCall(PetscLogEventBegin(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL)); 1612 1613 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary)); 1614 1615 /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */ 1616 if (binary) { 1617 parentviewer = viewer; 1618 PetscCall(PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer)); 1619 } 1620 1621 if (rank == 0) { 1622 GmshFile gmsh[1]; 1623 char line[PETSC_MAX_PATH_LEN]; 1624 PetscBool match; 1625 1626 PetscCall(PetscArrayzero(gmsh, 1)); 1627 gmsh->viewer = viewer; 1628 gmsh->binary = binary; 1629 1630 PetscCall(GmshMeshCreate(&mesh)); 1631 1632 /* Read mesh format */ 1633 PetscCall(GmshReadSection(gmsh, line)); 1634 PetscCall(GmshExpect(gmsh, "$MeshFormat", line)); 1635 PetscCall(GmshReadMeshFormat(gmsh)); 1636 PetscCall(GmshReadEndSection(gmsh, "$EndMeshFormat", line)); 1637 1638 /* OPTIONAL Read mesh version (Neper only) */ 1639 PetscCall(GmshReadSection(gmsh, line)); 1640 PetscCall(GmshMatch(gmsh, "$MeshVersion", line, &match)); 1641 if (match) { 1642 PetscCall(GmshExpect(gmsh, "$MeshVersion", line)); 1643 PetscCall(GmshReadMeshVersion(gmsh)); 1644 PetscCall(GmshReadEndSection(gmsh, "$EndMeshVersion", line)); 1645 /* Initial read for entity section */ 1646 PetscCall(GmshReadSection(gmsh, line)); 1647 } 1648 1649 /* OPTIONAL Read mesh domain (Neper only) */ 1650 PetscCall(GmshMatch(gmsh, "$Domain", line, &match)); 1651 if (match) { 1652 PetscCall(GmshExpect(gmsh, "$Domain", line)); 1653 PetscCall(GmshReadMeshDomain(gmsh)); 1654 PetscCall(GmshReadEndSection(gmsh, "$EndDomain", line)); 1655 /* Initial read for entity section */ 1656 PetscCall(GmshReadSection(gmsh, line)); 1657 } 1658 1659 /* OPTIONAL Read physical names */ 1660 PetscCall(GmshMatch(gmsh, "$PhysicalNames", line, &match)); 1661 if (match) { 1662 PetscCall(GmshExpect(gmsh, "$PhysicalNames", line)); 1663 PetscCall(GmshReadPhysicalNames(gmsh, mesh)); 1664 PetscCall(GmshReadEndSection(gmsh, "$EndPhysicalNames", line)); 1665 /* Initial read for entity section */ 1666 PetscCall(GmshReadSection(gmsh, line)); 1667 } 1668 1669 /* OPTIONAL Read entities */ 1670 if (gmsh->fileFormat >= 40) { 1671 PetscCall(GmshMatch(gmsh, "$Entities", line, &match)); 1672 if (match) { 1673 PetscCall(GmshExpect(gmsh, "$Entities", line)); 1674 PetscCall(GmshReadEntities(gmsh, mesh)); 1675 PetscCall(GmshReadEndSection(gmsh, "$EndEntities", line)); 1676 /* Initial read for nodes section */ 1677 PetscCall(GmshReadSection(gmsh, line)); 1678 } 1679 } 1680 1681 /* Read nodes */ 1682 PetscCall(GmshExpect(gmsh, "$Nodes", line)); 1683 PetscCall(GmshReadNodes(gmsh, mesh)); 1684 PetscCall(GmshReadEndSection(gmsh, "$EndNodes", line)); 1685 1686 /* Read elements */ 1687 PetscCall(GmshReadSection(gmsh, line)); 1688 PetscCall(GmshExpect(gmsh, "$Elements", line)); 1689 PetscCall(GmshReadElements(gmsh, mesh)); 1690 PetscCall(GmshReadEndSection(gmsh, "$EndElements", line)); 1691 1692 /* Read periodic section (OPTIONAL) */ 1693 if (periodic) { 1694 PetscCall(GmshReadSection(gmsh, line)); 1695 PetscCall(GmshMatch(gmsh, "$Periodic", line, &periodic)); 1696 } 1697 if (periodic) { 1698 PetscCall(GmshExpect(gmsh, "$Periodic", line)); 1699 PetscCall(GmshReadPeriodic(gmsh, mesh)); 1700 PetscCall(GmshReadEndSection(gmsh, "$EndPeriodic", line)); 1701 } 1702 1703 PetscCall(PetscFree(gmsh->wbuf)); 1704 PetscCall(PetscFree(gmsh->sbuf)); 1705 PetscCall(PetscFree(gmsh->nbuf)); 1706 1707 dim = mesh->dim; 1708 order = mesh->order; 1709 numNodes = mesh->numNodes; 1710 numElems = mesh->numElems; 1711 numVerts = mesh->numVerts; 1712 numCells = mesh->numCells; 1713 1714 { 1715 GmshElement *elemA = mesh->numCells > 0 ? mesh->elements : NULL; 1716 GmshElement *elemB = PetscSafePointerPlusOffset(elemA, mesh->numCells - 1); 1717 int ptA = elemA ? GmshCellMap[elemA->cellType].polytope : -1; 1718 int ptB = elemB ? GmshCellMap[elemB->cellType].polytope : -1; 1719 isSimplex = (ptA == GMSH_QUA || ptA == GMSH_HEX) ? PETSC_FALSE : PETSC_TRUE; 1720 isHybrid = (ptA == ptB) ? PETSC_FALSE : PETSC_TRUE; 1721 hasTetra = (ptA == GMSH_TET) ? PETSC_TRUE : PETSC_FALSE; 1722 } 1723 } 1724 1725 if (parentviewer) PetscCall(PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer)); 1726 1727 { 1728 int buf[6]; 1729 1730 buf[0] = (int)dim; 1731 buf[1] = (int)order; 1732 buf[2] = periodic; 1733 buf[3] = isSimplex; 1734 buf[4] = isHybrid; 1735 buf[5] = hasTetra; 1736 1737 PetscCallMPI(MPI_Bcast(buf, 6, MPI_INT, 0, comm)); 1738 1739 dim = buf[0]; 1740 order = buf[1]; 1741 periodic = buf[2] ? PETSC_TRUE : PETSC_FALSE; 1742 isSimplex = buf[3] ? PETSC_TRUE : PETSC_FALSE; 1743 isHybrid = buf[4] ? PETSC_TRUE : PETSC_FALSE; 1744 hasTetra = buf[5] ? PETSC_TRUE : PETSC_FALSE; 1745 } 1746 1747 if (!highOrderSet) highOrder = (order > 1) ? PETSC_TRUE : PETSC_FALSE; 1748 PetscCheck(!highOrder || !isHybrid, comm, PETSC_ERR_SUP, "No support for discretization on hybrid meshes yet"); 1749 1750 /* We do not want this label automatically computed, instead we fill it here */ 1751 PetscCall(DMCreateLabel(*dm, "celltype")); 1752 1753 /* Allocate the cell-vertex mesh */ 1754 PetscCall(DMPlexSetChart(*dm, 0, numCells + numVerts)); 1755 for (cell = 0; cell < numCells; ++cell) { 1756 GmshElement *elem = mesh->elements + cell; 1757 DMPolytopeType ctype = DMPolytopeTypeFromGmsh(elem->cellType); 1758 if (hybrid && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR; 1759 PetscCall(DMPlexSetConeSize(*dm, cell, elem->numVerts)); 1760 PetscCall(DMPlexSetCellType(*dm, cell, ctype)); 1761 } 1762 for (v = numCells; v < numCells + numVerts; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT)); 1763 PetscCall(DMSetUp(*dm)); 1764 1765 /* Add cell-vertex connections */ 1766 for (cell = 0; cell < numCells; ++cell) { 1767 GmshElement *elem = mesh->elements + cell; 1768 for (v = 0; v < elem->numVerts; ++v) { 1769 const PetscInt nn = elem->nodes[v]; 1770 const PetscInt vv = mesh->vertexMap[nn]; 1771 cone[v] = numCells + vv; 1772 } 1773 PetscCall(DMPlexReorderCell(*dm, cell, cone)); 1774 PetscCall(DMPlexSetCone(*dm, cell, cone)); 1775 } 1776 1777 PetscCall(DMSetDimension(*dm, dim)); 1778 PetscCall(DMPlexSymmetrize(*dm)); 1779 PetscCall(DMPlexStratify(*dm)); 1780 if (interpolate) { 1781 DM idm; 1782 1783 PetscCall(DMPlexInterpolate(*dm, &idm)); 1784 PetscCall(DMDestroy(dm)); 1785 *dm = idm; 1786 } 1787 PetscCall(DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd)); 1788 1789 if (rank == 0) { 1790 const PetscInt Nr = useregions ? mesh->numRegions : 0; 1791 1792 PetscCall(PetscCalloc1(Nr, ®ionSets)); 1793 for (cell = 0, e = 0; e < numElems; ++e) { 1794 GmshElement *elem = mesh->elements + e; 1795 1796 /* Create cell sets */ 1797 if (elem->dim == dim && dim > 0) { 1798 if (elem->numTags > 0) { 1799 PetscInt Nt = elem->numTags, t, r; 1800 1801 for (t = 0; t < Nt; ++t) { 1802 const PetscInt tag = elem->tags[t]; 1803 const PetscBool generic = usegeneric && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE; 1804 1805 if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", cell, tag)); 1806 for (r = 0; r < Nr; ++r) { 1807 if (mesh->regionDims[r] != dim) continue; 1808 if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], cell, tag)); 1809 } 1810 } 1811 } 1812 cell++; 1813 } 1814 1815 /* Create face sets */ 1816 if (elem->numTags && interpolate && elem->dim == dim - 1) { 1817 PetscInt joinSize; 1818 const PetscInt *join = NULL; 1819 PetscInt Nt = elem->numTags, pdepth; 1820 1821 /* Find the relevant facet with vertex joins */ 1822 for (v = 0; v < elem->numVerts; ++v) { 1823 const PetscInt nn = elem->nodes[v]; 1824 const PetscInt vv = mesh->vertexMap[nn]; 1825 cone[v] = vStart + vv; 1826 } 1827 PetscCall(DMPlexGetFullJoin(*dm, elem->numVerts, cone, &joinSize, &join)); 1828 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); 1829 PetscCall(DMPlexGetPointDepth(*dm, join[0], &pdepth)); 1830 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); 1831 for (PetscInt t = 0; t < Nt; ++t) { 1832 const PetscInt tag = elem->tags[t]; 1833 const PetscBool generic = usegeneric && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE; 1834 1835 if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], tag)); 1836 for (PetscInt r = 0; r < Nr; ++r) { 1837 if (mesh->regionDims[r] != dim - 1) continue; 1838 if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], join[0], tag)); 1839 } 1840 } 1841 PetscCall(DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join)); 1842 } 1843 1844 /* Create edge sets */ 1845 if (elem->numTags && interpolate && dim > 2 && elem->dim == 1) { 1846 PetscInt joinSize; 1847 const PetscInt *join = NULL; 1848 PetscInt Nt = elem->numTags, t, r; 1849 1850 /* Find the relevant edge with vertex joins */ 1851 for (v = 0; v < elem->numVerts; ++v) { 1852 const PetscInt nn = elem->nodes[v]; 1853 const PetscInt vv = mesh->vertexMap[nn]; 1854 cone[v] = vStart + vv; 1855 } 1856 PetscCall(DMPlexGetJoin(*dm, elem->numVerts, cone, &joinSize, &join)); 1857 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); 1858 for (t = 0; t < Nt; ++t) { 1859 const PetscInt tag = elem->tags[t]; 1860 const PetscBool generic = usegeneric && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE; 1861 1862 if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &edgeSets, "Edge Sets", join[0], tag)); 1863 for (r = 0; r < Nr; ++r) { 1864 if (mesh->regionDims[r] != 1) continue; 1865 if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], join[0], tag)); 1866 } 1867 } 1868 PetscCall(DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join)); 1869 } 1870 1871 /* Create vertex sets */ 1872 if (elem->numTags && elem->dim == 0 && markvertices) { 1873 const PetscInt nn = elem->nodes[0]; 1874 const PetscInt vv = mesh->vertexMap[nn]; 1875 const PetscInt tag = elem->tags[0]; 1876 PetscInt r; 1877 1878 if (vv < 0) continue; 1879 if (usegeneric) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag)); 1880 for (r = 0; r < Nr; ++r) { 1881 if (mesh->regionDims[r] != 0) continue; 1882 if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag)); 1883 } 1884 } 1885 } 1886 if (markvertices) { 1887 for (v = 0; v < numNodes; ++v) { 1888 const PetscInt vv = mesh->vertexMap[v]; 1889 const PetscInt *tags = &mesh->nodelist->tag[v * GMSH_MAX_TAGS]; 1890 PetscInt r, t; 1891 1892 if (vv < 0) continue; 1893 for (t = 0; t < GMSH_MAX_TAGS; ++t) { 1894 const PetscInt tag = tags[t]; 1895 const PetscBool generic = usegeneric && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE; 1896 1897 if (tag == -1) continue; 1898 if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag)); 1899 for (r = 0; r < Nr; ++r) { 1900 if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag)); 1901 } 1902 } 1903 } 1904 } 1905 PetscCall(PetscFree(regionSets)); 1906 } 1907 1908 { /* Create Cell/Face/Edge/Vertex Sets labels at all processes */ 1909 enum { 1910 n = 5 1911 }; 1912 PetscBool flag[n]; 1913 1914 flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE; 1915 flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE; 1916 flag[2] = edgeSets ? PETSC_TRUE : PETSC_FALSE; 1917 flag[3] = vertSets ? PETSC_TRUE : PETSC_FALSE; 1918 flag[4] = marker ? PETSC_TRUE : PETSC_FALSE; 1919 PetscCallMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm)); 1920 if (flag[0]) PetscCall(DMCreateLabel(*dm, "Cell Sets")); 1921 if (flag[1]) PetscCall(DMCreateLabel(*dm, "Face Sets")); 1922 if (flag[2]) PetscCall(DMCreateLabel(*dm, "Edge Sets")); 1923 if (flag[3]) PetscCall(DMCreateLabel(*dm, "Vertex Sets")); 1924 if (flag[4]) PetscCall(DMCreateLabel(*dm, "marker")); 1925 } 1926 1927 if (periodic) { 1928 PetscCall(PetscBTCreate(numVerts, &periodicVerts)); 1929 for (n = 0; n < numNodes; ++n) { 1930 if (mesh->vertexMap[n] >= 0) { 1931 if (PetscUnlikely(mesh->periodMap[n] != n)) { 1932 PetscInt m = mesh->periodMap[n]; 1933 PetscCall(PetscBTSet(periodicVerts, mesh->vertexMap[n])); 1934 PetscCall(PetscBTSet(periodicVerts, mesh->vertexMap[m])); 1935 } 1936 } 1937 } 1938 PetscCall(DMGetCoordinateDM(*dm, &cdm)); 1939 PetscCall(PetscMalloc1(maxHeight + 1, &periodicCells)); 1940 for (PetscInt h = 0; h <= maxHeight; ++h) { 1941 PetscInt pStart, pEnd; 1942 1943 PetscCall(DMPlexGetHeightStratum(*dm, h, &pStart, &pEnd)); 1944 PetscCall(PetscBTCreate(pEnd - pStart, &periodicCells[h])); 1945 for (PetscInt p = pStart; p < pEnd; ++p) { 1946 PetscInt *closure = NULL; 1947 PetscInt Ncl; 1948 1949 PetscCall(DMPlexGetTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure)); 1950 for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) { 1951 if (closure[cl] >= vStart && closure[cl] < vEnd) { 1952 if (PetscUnlikely(PetscBTLookup(periodicVerts, closure[cl] - vStart))) { 1953 PetscCall(PetscBTSet(periodicCells[h], p - pStart)); 1954 break; 1955 } 1956 } 1957 } 1958 PetscCall(DMPlexRestoreTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure)); 1959 } 1960 } 1961 } 1962 1963 /* Setup coordinate DM */ 1964 if (coordDim < 0) coordDim = dim; 1965 PetscCall(DMSetCoordinateDim(*dm, coordDim)); 1966 PetscCall(DMGetCoordinateDM(*dm, &cdm)); 1967 if (highOrder) { 1968 PetscFE fe; 1969 PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE; 1970 PetscDTNodeType nodeType = PETSCDTNODES_EQUISPACED; 1971 1972 if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */ 1973 1974 PetscCall(GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe)); 1975 PetscCall(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view")); 1976 PetscCall(DMSetField(cdm, 0, NULL, (PetscObject)fe)); 1977 PetscCall(PetscFEDestroy(&fe)); 1978 PetscCall(DMCreateDS(cdm)); 1979 } 1980 if (periodic) { 1981 PetscCall(DMClone(cdm, &cdmCell)); 1982 PetscCall(DMSetCellCoordinateDM(*dm, cdmCell)); 1983 } 1984 1985 /* Create coordinates */ 1986 if (highOrder) { 1987 PetscInt maxDof = GmshNumNodes_HEX(order) * coordDim; 1988 double *coords = mesh ? mesh->nodelist->xyz : NULL; 1989 PetscSection section; 1990 PetscScalar *cellCoords; 1991 1992 PetscCall(DMSetLocalSection(cdm, NULL)); 1993 PetscCall(DMGetLocalSection(cdm, &cs)); 1994 PetscCall(PetscSectionClone(cs, §ion)); 1995 PetscCall(DMPlexSetClosurePermutationTensor(cdm, 0, section)); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */ 1996 1997 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 1998 PetscCall(PetscMalloc1(maxDof, &cellCoords)); 1999 for (cell = 0; cell < numCells; ++cell) { 2000 GmshElement *elem = mesh->elements + cell; 2001 const int *lexorder = GmshCellMap[elem->cellType].lexorder(); 2002 int s = 0; 2003 for (n = 0; n < elem->numNodes; ++n) { 2004 while (lexorder[n + s] < 0) ++s; 2005 const PetscInt node = elem->nodes[lexorder[n + s]]; 2006 for (d = 0; d < coordDim; ++d) cellCoords[(n + s) * coordDim + d] = (PetscReal)coords[node * 3 + d]; 2007 } 2008 if (s) { 2009 /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */ 2010 PetscReal quaCenterWeights[9] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25}; 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 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}; 2013 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}; 2014 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}; 2015 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}; 2016 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}; 2017 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}; 2018 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}; 2019 PetscReal *sdWeights2[9] = {NULL, NULL, NULL, NULL, quaCenterWeights, NULL, NULL, NULL, NULL}; 2020 PetscReal *sdWeights3[27] = {NULL, NULL, NULL, NULL, hexBottomWeights, NULL, NULL, NULL, NULL, NULL, hexFrontWeights, NULL, hexLeftWeights, hexCenterWeights, hexRightWeights, NULL, hexBackWeights, NULL, 2021 NULL, NULL, NULL, NULL, hexTopWeights, NULL, NULL, NULL, NULL}; 2022 PetscReal **sdWeights[4] = {NULL, NULL, sdWeights2, sdWeights3}; 2023 2024 /* Missing entries in serendipity cell, only works for 8-node quad and 20-node hex */ 2025 for (n = 0; n < elem->numNodes + s; ++n) { 2026 if (lexorder[n] >= 0) continue; 2027 for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] = 0.0; 2028 for (int bn = 0; bn < elem->numNodes + s; ++bn) { 2029 if (lexorder[bn] < 0) continue; 2030 const PetscReal *weights = sdWeights[coordDim][n]; 2031 const PetscInt bnode = elem->nodes[lexorder[bn]]; 2032 for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] += weights[bn] * (PetscReal)coords[bnode * 3 + d]; 2033 } 2034 } 2035 } 2036 PetscCall(DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES)); 2037 } 2038 PetscCall(PetscSectionDestroy(§ion)); 2039 PetscCall(PetscFree(cellCoords)); 2040 } else { 2041 PetscInt *nodeMap; 2042 double *coords = mesh ? mesh->nodelist->xyz : NULL; 2043 PetscScalar *pointCoords; 2044 2045 PetscCall(DMGetCoordinateSection(*dm, &cs)); 2046 PetscCall(PetscSectionSetNumFields(cs, 1)); 2047 PetscCall(PetscSectionSetFieldComponents(cs, 0, coordDim)); 2048 PetscCall(PetscSectionSetChart(cs, numCells, numCells + numVerts)); 2049 for (v = numCells; v < numCells + numVerts; ++v) { 2050 PetscCall(PetscSectionSetDof(cs, v, coordDim)); 2051 PetscCall(PetscSectionSetFieldDof(cs, v, 0, coordDim)); 2052 } 2053 PetscCall(PetscSectionSetUp(cs)); 2054 2055 /* We need to localize coordinates on cells */ 2056 if (periodic) { 2057 PetscInt newStart = PETSC_INT_MAX, newEnd = -1, pStart, pEnd; 2058 2059 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)cdmCell), &csCell)); 2060 PetscCall(PetscSectionSetNumFields(csCell, 1)); 2061 PetscCall(PetscSectionSetFieldComponents(csCell, 0, coordDim)); 2062 for (PetscInt h = 0; h <= maxHeight; h++) { 2063 PetscCall(DMPlexGetHeightStratum(cdmCell, h, &pStart, &pEnd)); 2064 newStart = PetscMin(newStart, pStart); 2065 newEnd = PetscMax(newEnd, pEnd); 2066 } 2067 PetscCall(PetscSectionSetChart(csCell, newStart, newEnd)); 2068 for (PetscInt h = 0; h <= maxHeight; h++) { 2069 PetscCall(DMPlexGetHeightStratum(cdmCell, h, &pStart, &pEnd)); 2070 for (PetscInt p = pStart; p < pEnd; ++p) { 2071 PetscInt *closure = NULL; 2072 PetscInt Ncl, Nv = 0; 2073 2074 if (PetscUnlikely(PetscBTLookup(periodicCells[h], p - pStart))) { 2075 PetscCall(DMPlexGetTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure)); 2076 for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) { 2077 if (closure[cl] >= vStart && closure[cl] < vEnd) ++Nv; 2078 } 2079 PetscCall(DMPlexRestoreTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure)); 2080 PetscCall(PetscSectionSetDof(csCell, p, Nv * coordDim)); 2081 PetscCall(PetscSectionSetFieldDof(csCell, p, 0, Nv * coordDim)); 2082 } 2083 } 2084 } 2085 PetscCall(PetscSectionSetUp(csCell)); 2086 PetscCall(DMSetCellCoordinateSection(*dm, PETSC_DETERMINE, csCell)); 2087 } 2088 2089 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 2090 PetscCall(VecGetArray(coordinates, &pointCoords)); 2091 PetscCall(PetscMalloc1(numVerts, &nodeMap)); 2092 for (n = 0; n < numNodes; n++) 2093 if (mesh->vertexMap[n] >= 0) nodeMap[mesh->vertexMap[n]] = n; 2094 for (v = 0; v < numVerts; ++v) { 2095 PetscInt off, node = nodeMap[v]; 2096 2097 PetscCall(PetscSectionGetOffset(cs, numCells + v, &off)); 2098 for (d = 0; d < coordDim; ++d) pointCoords[off + d] = (PetscReal)coords[node * 3 + d]; 2099 } 2100 PetscCall(VecRestoreArray(coordinates, &pointCoords)); 2101 2102 if (periodic) { 2103 PetscInt cStart, cEnd; 2104 2105 PetscCall(DMPlexGetHeightStratum(cdmCell, 0, &cStart, &cEnd)); 2106 PetscCall(DMCreateLocalVector(cdmCell, &coordinatesCell)); 2107 PetscCall(VecGetArray(coordinatesCell, &pointCoords)); 2108 for (PetscInt c = cStart; c < cEnd; ++c) { 2109 GmshElement *elem = mesh->elements + c - cStart; 2110 PetscInt *closure = NULL; 2111 PetscInt verts[8]; 2112 PetscInt Ncl, Nv = 0; 2113 2114 for (PetscInt v = 0; v < elem->numVerts; ++v) cone[v] = elem->nodes[v]; 2115 PetscCall(DMPlexReorderCell(cdmCell, c, cone)); 2116 PetscCall(DMPlexGetTransitiveClosure(cdmCell, c, PETSC_TRUE, &Ncl, &closure)); 2117 for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) { 2118 if (closure[cl] >= vStart && closure[cl] < vEnd) verts[Nv++] = closure[cl]; 2119 } 2120 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); 2121 for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) { 2122 const PetscInt point = closure[cl]; 2123 PetscInt hStart, h; 2124 2125 PetscCall(DMPlexGetPointHeight(cdmCell, point, &h)); 2126 if (h > maxHeight) continue; 2127 PetscCall(DMPlexGetHeightStratum(cdmCell, h, &hStart, NULL)); 2128 if (PetscUnlikely(PetscBTLookup(periodicCells[h], point - hStart))) { 2129 PetscInt *pclosure = NULL; 2130 PetscInt Npcl, off, v; 2131 2132 PetscCall(PetscSectionGetOffset(csCell, point, &off)); 2133 PetscCall(DMPlexGetTransitiveClosure(cdmCell, point, PETSC_TRUE, &Npcl, &pclosure)); 2134 for (PetscInt pcl = 0; pcl < Npcl * 2; pcl += 2) { 2135 if (pclosure[pcl] >= vStart && pclosure[pcl] < vEnd) { 2136 for (v = 0; v < Nv; ++v) 2137 if (verts[v] == pclosure[pcl]) break; 2138 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); 2139 for (PetscInt d = 0; d < coordDim; ++d) pointCoords[off++] = (PetscReal)coords[cone[v] * 3 + d]; 2140 ++v; 2141 } 2142 } 2143 PetscCall(DMPlexRestoreTransitiveClosure(cdmCell, point, PETSC_TRUE, &Npcl, &pclosure)); 2144 } 2145 } 2146 PetscCall(DMPlexRestoreTransitiveClosure(cdmCell, c, PETSC_TRUE, &Ncl, &closure)); 2147 } 2148 PetscCall(VecSetBlockSize(coordinatesCell, coordDim)); 2149 PetscCall(VecRestoreArray(coordinatesCell, &pointCoords)); 2150 PetscCall(DMSetCellCoordinatesLocal(*dm, coordinatesCell)); 2151 PetscCall(VecDestroy(&coordinatesCell)); 2152 } 2153 PetscCall(PetscFree(nodeMap)); 2154 PetscCall(PetscSectionDestroy(&csCell)); 2155 PetscCall(DMDestroy(&cdmCell)); 2156 } 2157 2158 PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates")); 2159 PetscCall(VecSetBlockSize(coordinates, coordDim)); 2160 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 2161 PetscCall(VecDestroy(&coordinates)); 2162 2163 PetscCall(GmshMeshDestroy(&mesh)); 2164 PetscCall(PetscBTDestroy(&periodicVerts)); 2165 if (periodic) { 2166 for (PetscInt h = 0; h <= maxHeight; ++h) PetscCall(PetscBTDestroy(periodicCells + h)); 2167 PetscCall(PetscFree(periodicCells)); 2168 } 2169 2170 if (highOrder && project) { 2171 PetscFE fe; 2172 const char prefix[] = "dm_plex_gmsh_project_"; 2173 PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE; 2174 PetscDTNodeType nodeType = PETSCDTNODES_GAUSSJACOBI; 2175 2176 if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */ 2177 PetscCall(GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe)); 2178 PetscCall(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view")); 2179 PetscCall(DMSetCoordinateDisc(*dm, fe, PETSC_TRUE)); 2180 PetscCall(PetscFEDestroy(&fe)); 2181 } 2182 2183 PetscCall(PetscLogEventEnd(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL)); 2184 PetscFunctionReturn(PETSC_SUCCESS); 2185 } 2186