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