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 0; 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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", Section); 271 PetscFunctionReturn(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 367 } 368 369 static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities) 370 { 371 PetscInt dim; 372 373 PetscFunctionBegin; 374 if (!*entities) PetscFunctionReturn(0); 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(0); 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(0); 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(0); 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(0); 417 } 418 419 static PetscErrorCode GmshNodesDestroy(GmshNodes **nodes) 420 { 421 PetscFunctionBegin; 422 if (!*nodes) PetscFunctionReturn(0); 423 PetscCall(PetscFree((*nodes)->id)); 424 PetscCall(PetscFree((*nodes)->xyz)); 425 PetscCall(PetscFree((*nodes)->tag)); 426 PetscCall(PetscFree((*nodes))); 427 PetscFunctionReturn(0); 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(0); 446 } 447 448 static PetscErrorCode GmshElementsDestroy(GmshElement **elements) 449 { 450 PetscFunctionBegin; 451 if (!*elements) PetscFunctionReturn(0); 452 PetscCall(PetscFree(*elements)); 453 PetscFunctionReturn(0); 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(0); 481 } 482 483 static PetscErrorCode GmshMeshDestroy(GmshMesh **mesh) 484 { 485 PetscInt r; 486 487 PetscFunctionBegin; 488 if (!*mesh) PetscFunctionReturn(0); 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(0); 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(0); 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(0); 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(0); 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 unitilized 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 1083 } 1084 1085 /* 1086 PhysicalNames 1087 numPhysicalNames(ASCII int) 1088 dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max) 1089 ... 1090 $EndPhysicalNames 1091 */ 1092 static PetscErrorCode GmshReadPhysicalNames(GmshFile *gmsh, GmshMesh *mesh) 1093 { 1094 char line[PETSC_MAX_PATH_LEN], name[128 + 2], *p, *q, *r; 1095 int snum, region, dim, tag; 1096 1097 PetscFunctionBegin; 1098 PetscCall(GmshReadString(gmsh, line, 1)); 1099 snum = sscanf(line, "%d", ®ion); 1100 mesh->numRegions = region; 1101 PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 1102 PetscCall(PetscMalloc3(mesh->numRegions, &mesh->regionDims, mesh->numRegions, &mesh->regionTags, mesh->numRegions, &mesh->regionNames)); 1103 for (region = 0; region < mesh->numRegions; ++region) { 1104 PetscCall(GmshReadString(gmsh, line, 2)); 1105 snum = sscanf(line, "%d %d", &dim, &tag); 1106 PetscCheck(snum == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 1107 PetscCall(GmshReadString(gmsh, line, -(PetscInt)sizeof(line))); 1108 PetscCall(PetscStrchr(line, '"', &p)); 1109 PetscCheck(p, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 1110 PetscCall(PetscStrrchr(line, '"', &q)); 1111 PetscCheck(q != p, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 1112 PetscCall(PetscStrrchr(line, ':', &r)); 1113 if (p != r) q = r; 1114 PetscCall(PetscStrncpy(name, p + 1, (size_t)(q - p - 1))); 1115 mesh->regionDims[region] = dim; 1116 mesh->regionTags[region] = tag; 1117 PetscCall(PetscStrallocpy(name, &mesh->regionNames[region])); 1118 } 1119 PetscFunctionReturn(0); 1120 } 1121 1122 static PetscErrorCode GmshReadEntities(GmshFile *gmsh, GmshMesh *mesh) 1123 { 1124 PetscFunctionBegin; 1125 switch (gmsh->fileFormat) { 1126 case 41: 1127 PetscCall(GmshReadEntities_v41(gmsh, mesh)); 1128 break; 1129 default: 1130 PetscCall(GmshReadEntities_v40(gmsh, mesh)); 1131 break; 1132 } 1133 PetscFunctionReturn(0); 1134 } 1135 1136 static PetscErrorCode GmshReadNodes(GmshFile *gmsh, GmshMesh *mesh) 1137 { 1138 PetscFunctionBegin; 1139 switch (gmsh->fileFormat) { 1140 case 41: 1141 PetscCall(GmshReadNodes_v41(gmsh, mesh)); 1142 break; 1143 case 40: 1144 PetscCall(GmshReadNodes_v40(gmsh, mesh)); 1145 break; 1146 default: 1147 PetscCall(GmshReadNodes_v22(gmsh, mesh)); 1148 break; 1149 } 1150 1151 { /* Gmsh v2.2/v4.0 does not provide min/max node tags */ 1152 if (mesh->numNodes > 0 && gmsh->nodeEnd >= gmsh->nodeStart) { 1153 PetscInt tagMin = PETSC_MAX_INT, tagMax = PETSC_MIN_INT, n; 1154 GmshNodes *nodes = mesh->nodelist; 1155 for (n = 0; n < mesh->numNodes; ++n) { 1156 const PetscInt tag = nodes->id[n]; 1157 tagMin = PetscMin(tag, tagMin); 1158 tagMax = PetscMax(tag, tagMax); 1159 } 1160 gmsh->nodeStart = tagMin; 1161 gmsh->nodeEnd = tagMax + 1; 1162 } 1163 } 1164 1165 { /* Support for sparse node tags */ 1166 PetscInt n, t; 1167 GmshNodes *nodes = mesh->nodelist; 1168 PetscCall(PetscMalloc1(gmsh->nodeEnd - gmsh->nodeStart, &gmsh->nbuf)); 1169 for (t = 0; t < gmsh->nodeEnd - gmsh->nodeStart; ++t) gmsh->nbuf[t] = PETSC_MIN_INT; 1170 gmsh->nodeMap = gmsh->nbuf - gmsh->nodeStart; 1171 for (n = 0; n < mesh->numNodes; ++n) { 1172 const PetscInt tag = nodes->id[n]; 1173 PetscCheck(gmsh->nodeMap[tag] < 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Repeated node tag %" PetscInt_FMT, tag); 1174 gmsh->nodeMap[tag] = n; 1175 } 1176 } 1177 PetscFunctionReturn(0); 1178 } 1179 1180 static PetscErrorCode GmshReadElements(GmshFile *gmsh, GmshMesh *mesh) 1181 { 1182 PetscFunctionBegin; 1183 switch (gmsh->fileFormat) { 1184 case 41: 1185 PetscCall(GmshReadElements_v41(gmsh, mesh)); 1186 break; 1187 case 40: 1188 PetscCall(GmshReadElements_v40(gmsh, mesh)); 1189 break; 1190 default: 1191 PetscCall(GmshReadElements_v22(gmsh, mesh)); 1192 break; 1193 } 1194 1195 { /* Reorder elements by codimension and polytope type */ 1196 PetscInt ne = mesh->numElems; 1197 GmshElement *elements = mesh->elements; 1198 PetscInt keymap[GMSH_NUM_POLYTOPES], nk = 0; 1199 PetscInt offset[GMSH_NUM_POLYTOPES + 1], e, k; 1200 1201 for (k = 0; k < GMSH_NUM_POLYTOPES; ++k) keymap[k] = PETSC_MIN_INT; 1202 PetscCall(PetscMemzero(offset, sizeof(offset))); 1203 1204 keymap[GMSH_TET] = nk++; 1205 keymap[GMSH_HEX] = nk++; 1206 keymap[GMSH_PRI] = nk++; 1207 keymap[GMSH_PYR] = nk++; 1208 keymap[GMSH_TRI] = nk++; 1209 keymap[GMSH_QUA] = nk++; 1210 keymap[GMSH_SEG] = nk++; 1211 keymap[GMSH_VTX] = nk++; 1212 1213 PetscCall(GmshElementsCreate(mesh->numElems, &mesh->elements)); 1214 #define key(eid) keymap[GmshCellMap[elements[(eid)].cellType].polytope] 1215 for (e = 0; e < ne; ++e) offset[1 + key(e)]++; 1216 for (k = 1; k < nk; ++k) offset[k] += offset[k - 1]; 1217 for (e = 0; e < ne; ++e) mesh->elements[offset[key(e)]++] = elements[e]; 1218 #undef key 1219 PetscCall(GmshElementsDestroy(&elements)); 1220 } 1221 1222 { /* Mesh dimension and order */ 1223 GmshElement *elem = mesh->numElems ? mesh->elements : NULL; 1224 mesh->dim = elem ? GmshCellMap[elem->cellType].dim : 0; 1225 mesh->order = elem ? GmshCellMap[elem->cellType].order : 0; 1226 } 1227 1228 { 1229 PetscBT vtx; 1230 PetscInt dim = mesh->dim, e, n, v; 1231 1232 PetscCall(PetscBTCreate(mesh->numNodes, &vtx)); 1233 1234 /* Compute number of cells and set of vertices */ 1235 mesh->numCells = 0; 1236 for (e = 0; e < mesh->numElems; ++e) { 1237 GmshElement *elem = mesh->elements + e; 1238 if (elem->dim == dim && dim > 0) mesh->numCells++; 1239 for (v = 0; v < elem->numVerts; v++) PetscCall(PetscBTSet(vtx, elem->nodes[v])); 1240 } 1241 1242 /* Compute numbering for vertices */ 1243 mesh->numVerts = 0; 1244 PetscCall(PetscMalloc1(mesh->numNodes, &mesh->vertexMap)); 1245 for (n = 0; n < mesh->numNodes; ++n) mesh->vertexMap[n] = PetscBTLookup(vtx, n) ? mesh->numVerts++ : PETSC_MIN_INT; 1246 1247 PetscCall(PetscBTDestroy(&vtx)); 1248 } 1249 PetscFunctionReturn(0); 1250 } 1251 1252 static PetscErrorCode GmshReadPeriodic(GmshFile *gmsh, GmshMesh *mesh) 1253 { 1254 PetscInt n; 1255 1256 PetscFunctionBegin; 1257 PetscCall(PetscMalloc1(mesh->numNodes, &mesh->periodMap)); 1258 for (n = 0; n < mesh->numNodes; ++n) mesh->periodMap[n] = n; 1259 switch (gmsh->fileFormat) { 1260 case 41: 1261 PetscCall(GmshReadPeriodic_v41(gmsh, mesh->periodMap)); 1262 break; 1263 default: 1264 PetscCall(GmshReadPeriodic_v40(gmsh, mesh->periodMap)); 1265 break; 1266 } 1267 1268 /* Find canonical primary nodes */ 1269 for (n = 0; n < mesh->numNodes; ++n) 1270 while (mesh->periodMap[n] != mesh->periodMap[mesh->periodMap[n]]) mesh->periodMap[n] = mesh->periodMap[mesh->periodMap[n]]; 1271 1272 /* Renumber vertices (filter out correspondings) */ 1273 mesh->numVerts = 0; 1274 for (n = 0; n < mesh->numNodes; ++n) 1275 if (mesh->vertexMap[n] >= 0) /* is vertex */ 1276 if (mesh->periodMap[n] == n) /* is primary */ 1277 mesh->vertexMap[n] = mesh->numVerts++; 1278 for (n = 0; n < mesh->numNodes; ++n) 1279 if (mesh->vertexMap[n] >= 0) /* is vertex */ 1280 if (mesh->periodMap[n] != n) /* is corresponding */ 1281 mesh->vertexMap[n] = mesh->vertexMap[mesh->periodMap[n]]; 1282 PetscFunctionReturn(0); 1283 } 1284 1285 #define DM_POLYTOPE_VERTEX DM_POLYTOPE_POINT 1286 static const DMPolytopeType DMPolytopeMap[] = { 1287 /* GMSH_VTX */ DM_POLYTOPE_VERTEX, 1288 /* GMSH_SEG */ DM_POLYTOPE_SEGMENT, 1289 /* GMSH_TRI */ DM_POLYTOPE_TRIANGLE, 1290 /* GMSH_QUA */ DM_POLYTOPE_QUADRILATERAL, 1291 /* GMSH_TET */ DM_POLYTOPE_TETRAHEDRON, 1292 /* GMSH_HEX */ DM_POLYTOPE_HEXAHEDRON, 1293 /* GMSH_PRI */ DM_POLYTOPE_TRI_PRISM, 1294 /* GMSH_PYR */ DM_POLYTOPE_PYRAMID, DM_POLYTOPE_UNKNOWN}; 1295 1296 static inline DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt cellType) 1297 { 1298 return DMPolytopeMap[GmshCellMap[cellType].polytope]; 1299 } 1300 1301 static PetscErrorCode GmshCreateFE(MPI_Comm comm, const char prefix[], PetscBool isSimplex, PetscBool continuity, PetscDTNodeType nodeType, PetscInt dim, PetscInt Nc, PetscInt k, PetscFE *fem) 1302 { 1303 DM K; 1304 PetscSpace P; 1305 PetscDualSpace Q; 1306 PetscQuadrature q, fq; 1307 PetscBool isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1308 PetscBool endpoint = PETSC_TRUE; 1309 char name[32]; 1310 1311 PetscFunctionBegin; 1312 /* Create space */ 1313 PetscCall(PetscSpaceCreate(comm, &P)); 1314 PetscCall(PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL)); 1315 PetscCall(PetscSpacePolynomialSetTensor(P, isTensor)); 1316 PetscCall(PetscSpaceSetNumComponents(P, Nc)); 1317 PetscCall(PetscSpaceSetNumVariables(P, dim)); 1318 PetscCall(PetscSpaceSetDegree(P, k, PETSC_DETERMINE)); 1319 if (prefix) { 1320 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)P, prefix)); 1321 PetscCall(PetscSpaceSetFromOptions(P)); 1322 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)P, NULL)); 1323 PetscCall(PetscSpaceGetDegree(P, &k, NULL)); 1324 } 1325 PetscCall(PetscSpaceSetUp(P)); 1326 /* Create dual space */ 1327 PetscCall(PetscDualSpaceCreate(comm, &Q)); 1328 PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE)); 1329 PetscCall(PetscDualSpaceLagrangeSetTensor(Q, isTensor)); 1330 PetscCall(PetscDualSpaceLagrangeSetContinuity(Q, continuity)); 1331 PetscCall(PetscDualSpaceLagrangeSetNodeType(Q, nodeType, endpoint, 0)); 1332 PetscCall(PetscDualSpaceSetNumComponents(Q, Nc)); 1333 PetscCall(PetscDualSpaceSetOrder(Q, k)); 1334 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K)); 1335 PetscCall(PetscDualSpaceSetDM(Q, K)); 1336 PetscCall(DMDestroy(&K)); 1337 if (prefix) { 1338 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)Q, prefix)); 1339 PetscCall(PetscDualSpaceSetFromOptions(Q)); 1340 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)Q, NULL)); 1341 } 1342 PetscCall(PetscDualSpaceSetUp(Q)); 1343 /* Create quadrature */ 1344 if (isSimplex) { 1345 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, k + 1, -1, +1, &q)); 1346 PetscCall(PetscDTStroudConicalQuadrature(dim - 1, 1, k + 1, -1, +1, &fq)); 1347 } else { 1348 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, k + 1, -1, +1, &q)); 1349 PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, k + 1, -1, +1, &fq)); 1350 } 1351 /* Create finite element */ 1352 PetscCall(PetscFECreate(comm, fem)); 1353 PetscCall(PetscFESetType(*fem, PETSCFEBASIC)); 1354 PetscCall(PetscFESetNumComponents(*fem, Nc)); 1355 PetscCall(PetscFESetBasisSpace(*fem, P)); 1356 PetscCall(PetscFESetDualSpace(*fem, Q)); 1357 PetscCall(PetscFESetQuadrature(*fem, q)); 1358 PetscCall(PetscFESetFaceQuadrature(*fem, fq)); 1359 if (prefix) { 1360 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*fem, prefix)); 1361 PetscCall(PetscFESetFromOptions(*fem)); 1362 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*fem, NULL)); 1363 } 1364 PetscCall(PetscFESetUp(*fem)); 1365 /* Cleanup */ 1366 PetscCall(PetscSpaceDestroy(&P)); 1367 PetscCall(PetscDualSpaceDestroy(&Q)); 1368 PetscCall(PetscQuadratureDestroy(&q)); 1369 PetscCall(PetscQuadratureDestroy(&fq)); 1370 /* Set finite element name */ 1371 PetscCall(PetscSNPrintf(name, sizeof(name), "%s%" PetscInt_FMT, isSimplex ? "P" : "Q", k)); 1372 PetscCall(PetscFESetName(*fem, name)); 1373 PetscFunctionReturn(0); 1374 } 1375 1376 /*@C 1377 DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file 1378 1379 + comm - The MPI communicator 1380 . filename - Name of the Gmsh file 1381 - interpolate - Create faces and edges in the mesh 1382 1383 Output Parameter: 1384 . dm - The DM object representing the mesh 1385 1386 Level: beginner 1387 1388 .seealso: `DMPlexCreateFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()` 1389 @*/ 1390 PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 1391 { 1392 PetscViewer viewer; 1393 PetscMPIInt rank; 1394 int fileType; 1395 PetscViewerType vtype; 1396 1397 PetscFunctionBegin; 1398 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1399 1400 /* Determine Gmsh file type (ASCII or binary) from file header */ 1401 if (rank == 0) { 1402 GmshFile gmsh[1]; 1403 char line[PETSC_MAX_PATH_LEN]; 1404 int snum; 1405 float version; 1406 int fileFormat; 1407 1408 PetscCall(PetscArrayzero(gmsh, 1)); 1409 PetscCall(PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer)); 1410 PetscCall(PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII)); 1411 PetscCall(PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ)); 1412 PetscCall(PetscViewerFileSetName(gmsh->viewer, filename)); 1413 /* Read only the first two lines of the Gmsh file */ 1414 PetscCall(GmshReadSection(gmsh, line)); 1415 PetscCall(GmshExpect(gmsh, "$MeshFormat", line)); 1416 PetscCall(GmshReadString(gmsh, line, 2)); 1417 snum = sscanf(line, "%f %d", &version, &fileType); 1418 fileFormat = (int)roundf(version * 10); 1419 PetscCheck(snum == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line); 1420 PetscCheck(fileFormat >= 22, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version); 1421 PetscCheck((int)version != 3, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version); 1422 PetscCheck(fileFormat <= 41, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version); 1423 PetscCall(PetscViewerDestroy(&gmsh->viewer)); 1424 } 1425 PetscCallMPI(MPI_Bcast(&fileType, 1, MPI_INT, 0, comm)); 1426 vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY; 1427 1428 /* Create appropriate viewer and build plex */ 1429 PetscCall(PetscViewerCreate(comm, &viewer)); 1430 PetscCall(PetscViewerSetType(viewer, vtype)); 1431 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 1432 PetscCall(PetscViewerFileSetName(viewer, filename)); 1433 PetscCall(DMPlexCreateGmsh(comm, viewer, interpolate, dm)); 1434 PetscCall(PetscViewerDestroy(&viewer)); 1435 PetscFunctionReturn(0); 1436 } 1437 1438 /*@ 1439 DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer 1440 1441 Collective 1442 1443 Input Parameters: 1444 + comm - The MPI communicator 1445 . viewer - The Viewer associated with a Gmsh file 1446 - interpolate - Create faces and edges in the mesh 1447 1448 Output Parameter: 1449 . dm - The DM object representing the mesh 1450 1451 Options Database: 1452 + -dm_plex_gmsh_hybrid - Force triangular prisms to use tensor order 1453 . -dm_plex_gmsh_periodic - Read Gmsh periodic section and construct a periodic Plex 1454 . -dm_plex_gmsh_highorder - Generate high-order coordinates 1455 . -dm_plex_gmsh_project - Project high-order coordinates to a different space, use the prefix dm_plex_gmsh_project_ to define the space 1456 . -dm_plex_gmsh_use_regions - Generate labels with region names 1457 . -dm_plex_gmsh_mark_vertices - Add vertices to generated labels 1458 . -dm_plex_gmsh_multiple_tags - Allow multiple tags for default labels 1459 - -dm_plex_gmsh_spacedim <d> - Embedding space dimension, if different from topological dimension 1460 1461 Note: 1462 The Gmsh file format is described in http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format 1463 1464 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. 1465 1466 Level: beginner 1467 1468 .seealso: `DMPLEX`, `DMCreate()` 1469 @*/ 1470 PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm) 1471 { 1472 GmshMesh *mesh = NULL; 1473 PetscViewer parentviewer = NULL; 1474 PetscBT periodicVerts = NULL; 1475 PetscBT periodicCells = NULL; 1476 DM cdm, cdmCell = NULL; 1477 PetscSection cs, csCell = NULL; 1478 Vec coordinates, coordinatesCell; 1479 DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL, marker = NULL, *regionSets; 1480 PetscInt dim = 0, coordDim = -1, order = 0; 1481 PetscInt numNodes = 0, numElems = 0, numVerts = 0, numCells = 0; 1482 PetscInt cell, cone[8], e, n, v, d; 1483 PetscBool binary, useregions = PETSC_FALSE, markvertices = PETSC_FALSE, multipleTags = PETSC_FALSE; 1484 PetscBool hybrid = interpolate, periodic = PETSC_TRUE; 1485 PetscBool highOrder = PETSC_TRUE, highOrderSet, project = PETSC_FALSE; 1486 PetscBool isSimplex = PETSC_FALSE, isHybrid = PETSC_FALSE, hasTetra = PETSC_FALSE; 1487 PetscMPIInt rank; 1488 1489 PetscFunctionBegin; 1490 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1491 PetscObjectOptionsBegin((PetscObject)viewer); 1492 PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Gmsh options"); 1493 PetscCall(PetscOptionsBool("-dm_plex_gmsh_hybrid", "Force triangular prisms to use tensor order", "DMPlexCreateGmsh", hybrid, &hybrid, NULL)); 1494 PetscCall(PetscOptionsBool("-dm_plex_gmsh_periodic", "Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL)); 1495 PetscCall(PetscOptionsBool("-dm_plex_gmsh_highorder", "Generate high-order coordinates", "DMPlexCreateGmsh", highOrder, &highOrder, &highOrderSet)); 1496 PetscCall(PetscOptionsBool("-dm_plex_gmsh_project", "Project high-order coordinates to a different space", "DMPlexCreateGmsh", project, &project, NULL)); 1497 PetscCall(PetscOptionsBool("-dm_plex_gmsh_use_regions", "Generate labels with region names", "DMPlexCreateGmsh", useregions, &useregions, NULL)); 1498 PetscCall(PetscOptionsBool("-dm_plex_gmsh_mark_vertices", "Add vertices to generated labels", "DMPlexCreateGmsh", markvertices, &markvertices, NULL)); 1499 PetscCall(PetscOptionsBool("-dm_plex_gmsh_multiple_tags", "Allow multiple tags for default labels", "DMPlexCreateGmsh", multipleTags, &multipleTags, NULL)); 1500 PetscCall(PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE)); 1501 PetscOptionsHeadEnd(); 1502 PetscOptionsEnd(); 1503 1504 PetscCall(GmshCellInfoSetUp()); 1505 1506 PetscCall(DMCreate(comm, dm)); 1507 PetscCall(DMSetType(*dm, DMPLEX)); 1508 PetscCall(PetscLogEventBegin(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL)); 1509 1510 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary)); 1511 1512 /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */ 1513 if (binary) { 1514 parentviewer = viewer; 1515 PetscCall(PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer)); 1516 } 1517 1518 if (rank == 0) { 1519 GmshFile gmsh[1]; 1520 char line[PETSC_MAX_PATH_LEN]; 1521 PetscBool match; 1522 1523 PetscCall(PetscArrayzero(gmsh, 1)); 1524 gmsh->viewer = viewer; 1525 gmsh->binary = binary; 1526 1527 PetscCall(GmshMeshCreate(&mesh)); 1528 1529 /* Read mesh format */ 1530 PetscCall(GmshReadSection(gmsh, line)); 1531 PetscCall(GmshExpect(gmsh, "$MeshFormat", line)); 1532 PetscCall(GmshReadMeshFormat(gmsh)); 1533 PetscCall(GmshReadEndSection(gmsh, "$EndMeshFormat", line)); 1534 1535 /* OPTIONAL Read physical names */ 1536 PetscCall(GmshReadSection(gmsh, line)); 1537 PetscCall(GmshMatch(gmsh, "$PhysicalNames", line, &match)); 1538 if (match) { 1539 PetscCall(GmshExpect(gmsh, "$PhysicalNames", line)); 1540 PetscCall(GmshReadPhysicalNames(gmsh, mesh)); 1541 PetscCall(GmshReadEndSection(gmsh, "$EndPhysicalNames", line)); 1542 /* Initial read for entity section */ 1543 PetscCall(GmshReadSection(gmsh, line)); 1544 } 1545 1546 /* Read entities */ 1547 if (gmsh->fileFormat >= 40) { 1548 PetscCall(GmshExpect(gmsh, "$Entities", line)); 1549 PetscCall(GmshReadEntities(gmsh, mesh)); 1550 PetscCall(GmshReadEndSection(gmsh, "$EndEntities", line)); 1551 /* Initial read for nodes section */ 1552 PetscCall(GmshReadSection(gmsh, line)); 1553 } 1554 1555 /* Read nodes */ 1556 PetscCall(GmshExpect(gmsh, "$Nodes", line)); 1557 PetscCall(GmshReadNodes(gmsh, mesh)); 1558 PetscCall(GmshReadEndSection(gmsh, "$EndNodes", line)); 1559 1560 /* Read elements */ 1561 PetscCall(GmshReadSection(gmsh, line)); 1562 PetscCall(GmshExpect(gmsh, "$Elements", line)); 1563 PetscCall(GmshReadElements(gmsh, mesh)); 1564 PetscCall(GmshReadEndSection(gmsh, "$EndElements", line)); 1565 1566 /* Read periodic section (OPTIONAL) */ 1567 if (periodic) { 1568 PetscCall(GmshReadSection(gmsh, line)); 1569 PetscCall(GmshMatch(gmsh, "$Periodic", line, &periodic)); 1570 } 1571 if (periodic) { 1572 PetscCall(GmshExpect(gmsh, "$Periodic", line)); 1573 PetscCall(GmshReadPeriodic(gmsh, mesh)); 1574 PetscCall(GmshReadEndSection(gmsh, "$EndPeriodic", line)); 1575 } 1576 1577 PetscCall(PetscFree(gmsh->wbuf)); 1578 PetscCall(PetscFree(gmsh->sbuf)); 1579 PetscCall(PetscFree(gmsh->nbuf)); 1580 1581 dim = mesh->dim; 1582 order = mesh->order; 1583 numNodes = mesh->numNodes; 1584 numElems = mesh->numElems; 1585 numVerts = mesh->numVerts; 1586 numCells = mesh->numCells; 1587 1588 { 1589 GmshElement *elemA = mesh->numCells > 0 ? mesh->elements : NULL; 1590 GmshElement *elemB = elemA ? elemA + mesh->numCells - 1 : NULL; 1591 int ptA = elemA ? GmshCellMap[elemA->cellType].polytope : -1; 1592 int ptB = elemB ? GmshCellMap[elemB->cellType].polytope : -1; 1593 isSimplex = (ptA == GMSH_QUA || ptA == GMSH_HEX) ? PETSC_FALSE : PETSC_TRUE; 1594 isHybrid = (ptA == ptB) ? PETSC_FALSE : PETSC_TRUE; 1595 hasTetra = (ptA == GMSH_TET) ? PETSC_TRUE : PETSC_FALSE; 1596 } 1597 } 1598 1599 if (parentviewer) PetscCall(PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer)); 1600 1601 { 1602 int buf[6]; 1603 1604 buf[0] = (int)dim; 1605 buf[1] = (int)order; 1606 buf[2] = periodic; 1607 buf[3] = isSimplex; 1608 buf[4] = isHybrid; 1609 buf[5] = hasTetra; 1610 1611 PetscCallMPI(MPI_Bcast(buf, 6, MPI_INT, 0, comm)); 1612 1613 dim = buf[0]; 1614 order = buf[1]; 1615 periodic = buf[2] ? PETSC_TRUE : PETSC_FALSE; 1616 isSimplex = buf[3] ? PETSC_TRUE : PETSC_FALSE; 1617 isHybrid = buf[4] ? PETSC_TRUE : PETSC_FALSE; 1618 hasTetra = buf[5] ? PETSC_TRUE : PETSC_FALSE; 1619 } 1620 1621 if (!highOrderSet) highOrder = (order > 1) ? PETSC_TRUE : PETSC_FALSE; 1622 PetscCheck(!highOrder || !isHybrid, comm, PETSC_ERR_SUP, "No support for discretization on hybrid meshes yet"); 1623 1624 /* We do not want this label automatically computed, instead we fill it here */ 1625 PetscCall(DMCreateLabel(*dm, "celltype")); 1626 1627 /* Allocate the cell-vertex mesh */ 1628 PetscCall(DMPlexSetChart(*dm, 0, numCells + numVerts)); 1629 for (cell = 0; cell < numCells; ++cell) { 1630 GmshElement *elem = mesh->elements + cell; 1631 DMPolytopeType ctype = DMPolytopeTypeFromGmsh(elem->cellType); 1632 if (hybrid && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR; 1633 PetscCall(DMPlexSetConeSize(*dm, cell, elem->numVerts)); 1634 PetscCall(DMPlexSetCellType(*dm, cell, ctype)); 1635 } 1636 for (v = numCells; v < numCells + numVerts; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT)); 1637 PetscCall(DMSetUp(*dm)); 1638 1639 /* Add cell-vertex connections */ 1640 for (cell = 0; cell < numCells; ++cell) { 1641 GmshElement *elem = mesh->elements + cell; 1642 for (v = 0; v < elem->numVerts; ++v) { 1643 const PetscInt nn = elem->nodes[v]; 1644 const PetscInt vv = mesh->vertexMap[nn]; 1645 cone[v] = numCells + vv; 1646 } 1647 PetscCall(DMPlexReorderCell(*dm, cell, cone)); 1648 PetscCall(DMPlexSetCone(*dm, cell, cone)); 1649 } 1650 1651 PetscCall(DMSetDimension(*dm, dim)); 1652 PetscCall(DMPlexSymmetrize(*dm)); 1653 PetscCall(DMPlexStratify(*dm)); 1654 if (interpolate) { 1655 DM idm; 1656 1657 PetscCall(DMPlexInterpolate(*dm, &idm)); 1658 PetscCall(DMDestroy(dm)); 1659 *dm = idm; 1660 } 1661 1662 if (rank == 0) { 1663 const PetscInt Nr = useregions ? mesh->numRegions : 0; 1664 PetscInt vStart, vEnd; 1665 1666 PetscCall(PetscCalloc1(Nr, ®ionSets)); 1667 PetscCall(DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd)); 1668 for (cell = 0, e = 0; e < numElems; ++e) { 1669 GmshElement *elem = mesh->elements + e; 1670 1671 /* Create cell sets */ 1672 if (elem->dim == dim && dim > 0) { 1673 if (elem->numTags > 0) { 1674 PetscInt Nt = elem->numTags, t, r; 1675 1676 for (t = 0; t < Nt; ++t) { 1677 const PetscInt tag = elem->tags[t]; 1678 const PetscBool generic = !Nr && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE; 1679 1680 if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", cell, tag)); 1681 for (r = 0; r < Nr; ++r) { 1682 if (mesh->regionDims[r] != dim) continue; 1683 if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], cell, tag)); 1684 } 1685 } 1686 } 1687 cell++; 1688 } 1689 1690 /* Create face sets */ 1691 if (interpolate && elem->dim == dim - 1) { 1692 PetscInt joinSize; 1693 const PetscInt *join = NULL; 1694 PetscInt Nt = elem->numTags, t, r; 1695 1696 /* Find the relevant facet with vertex joins */ 1697 for (v = 0; v < elem->numVerts; ++v) { 1698 const PetscInt nn = elem->nodes[v]; 1699 const PetscInt vv = mesh->vertexMap[nn]; 1700 cone[v] = vStart + vv; 1701 } 1702 PetscCall(DMPlexGetFullJoin(*dm, elem->numVerts, cone, &joinSize, &join)); 1703 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); 1704 for (t = 0; t < Nt; ++t) { 1705 const PetscInt tag = elem->tags[t]; 1706 const PetscBool generic = !Nr && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE; 1707 1708 if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], tag)); 1709 for (r = 0; r < Nr; ++r) { 1710 if (mesh->regionDims[r] != dim - 1) continue; 1711 if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], join[0], tag)); 1712 } 1713 } 1714 PetscCall(DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join)); 1715 } 1716 1717 /* Create vertex sets */ 1718 if (elem->dim == 0 && markvertices) { 1719 if (elem->numTags > 0) { 1720 const PetscInt nn = elem->nodes[0]; 1721 const PetscInt vv = mesh->vertexMap[nn]; 1722 const PetscInt tag = elem->tags[0]; 1723 PetscInt r; 1724 1725 if (!Nr) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag)); 1726 for (r = 0; r < Nr; ++r) { 1727 if (mesh->regionDims[r] != 0) continue; 1728 if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag)); 1729 } 1730 } 1731 } 1732 } 1733 if (markvertices) { 1734 for (v = 0; v < numNodes; ++v) { 1735 const PetscInt vv = mesh->vertexMap[v]; 1736 const PetscInt *tags = &mesh->nodelist->tag[v * GMSH_MAX_TAGS]; 1737 PetscInt r, t; 1738 1739 for (t = 0; t < GMSH_MAX_TAGS; ++t) { 1740 const PetscInt tag = tags[t]; 1741 const PetscBool generic = !Nr && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE; 1742 1743 if (tag == -1) continue; 1744 if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag)); 1745 for (r = 0; r < Nr; ++r) { 1746 if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag)); 1747 } 1748 } 1749 } 1750 } 1751 PetscCall(PetscFree(regionSets)); 1752 } 1753 1754 { /* Create Cell/Face/Vertex Sets labels at all processes */ 1755 enum { 1756 n = 4 1757 }; 1758 PetscBool flag[n]; 1759 1760 flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE; 1761 flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE; 1762 flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE; 1763 flag[3] = marker ? PETSC_TRUE : PETSC_FALSE; 1764 PetscCallMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm)); 1765 if (flag[0]) PetscCall(DMCreateLabel(*dm, "Cell Sets")); 1766 if (flag[1]) PetscCall(DMCreateLabel(*dm, "Face Sets")); 1767 if (flag[2]) PetscCall(DMCreateLabel(*dm, "Vertex Sets")); 1768 if (flag[3]) PetscCall(DMCreateLabel(*dm, "marker")); 1769 } 1770 1771 if (periodic) { 1772 PetscCall(PetscBTCreate(numVerts, &periodicVerts)); 1773 for (n = 0; n < numNodes; ++n) { 1774 if (mesh->vertexMap[n] >= 0) { 1775 if (PetscUnlikely(mesh->periodMap[n] != n)) { 1776 PetscInt m = mesh->periodMap[n]; 1777 PetscCall(PetscBTSet(periodicVerts, mesh->vertexMap[n])); 1778 PetscCall(PetscBTSet(periodicVerts, mesh->vertexMap[m])); 1779 } 1780 } 1781 } 1782 PetscCall(PetscBTCreate(numCells, &periodicCells)); 1783 for (cell = 0; cell < numCells; ++cell) { 1784 GmshElement *elem = mesh->elements + cell; 1785 for (v = 0; v < elem->numVerts; ++v) { 1786 PetscInt nn = elem->nodes[v]; 1787 PetscInt vv = mesh->vertexMap[nn]; 1788 if (PetscUnlikely(PetscBTLookup(periodicVerts, vv))) { 1789 PetscCall(PetscBTSet(periodicCells, cell)); 1790 break; 1791 } 1792 } 1793 } 1794 } 1795 1796 /* Setup coordinate DM */ 1797 if (coordDim < 0) coordDim = dim; 1798 PetscCall(DMSetCoordinateDim(*dm, coordDim)); 1799 PetscCall(DMGetCoordinateDM(*dm, &cdm)); 1800 if (highOrder) { 1801 PetscFE fe; 1802 PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE; 1803 PetscDTNodeType nodeType = PETSCDTNODES_EQUISPACED; 1804 1805 if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */ 1806 1807 PetscCall(GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe)); 1808 PetscCall(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view")); 1809 PetscCall(DMSetField(cdm, 0, NULL, (PetscObject)fe)); 1810 PetscCall(PetscFEDestroy(&fe)); 1811 PetscCall(DMCreateDS(cdm)); 1812 } 1813 if (periodic) { 1814 PetscCall(DMClone(cdm, &cdmCell)); 1815 PetscCall(DMSetCellCoordinateDM(*dm, cdmCell)); 1816 } 1817 1818 /* Create coordinates */ 1819 if (highOrder) { 1820 PetscInt maxDof = GmshNumNodes_HEX(order) * coordDim; 1821 double *coords = mesh ? mesh->nodelist->xyz : NULL; 1822 PetscSection section; 1823 PetscScalar *cellCoords; 1824 1825 PetscCall(DMSetLocalSection(cdm, NULL)); 1826 PetscCall(DMGetLocalSection(cdm, &cs)); 1827 PetscCall(PetscSectionClone(cs, §ion)); 1828 PetscCall(DMPlexSetClosurePermutationTensor(cdm, 0, section)); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */ 1829 1830 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 1831 PetscCall(PetscMalloc1(maxDof, &cellCoords)); 1832 for (cell = 0; cell < numCells; ++cell) { 1833 GmshElement *elem = mesh->elements + cell; 1834 const int *lexorder = GmshCellMap[elem->cellType].lexorder(); 1835 int s = 0; 1836 for (n = 0; n < elem->numNodes; ++n) { 1837 while (lexorder[n + s] < 0) ++s; 1838 const PetscInt node = elem->nodes[lexorder[n + s]]; 1839 for (d = 0; d < coordDim; ++d) cellCoords[(n + s) * coordDim + d] = (PetscReal)coords[node * 3 + d]; 1840 } 1841 if (s) { 1842 /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */ 1843 PetscReal quaCenterWeights[9] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25}; 1844 /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */ 1845 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}; 1846 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}; 1847 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}; 1848 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}; 1849 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}; 1850 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}; 1851 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}; 1852 PetscReal *sdWeights2[9] = {NULL, NULL, NULL, NULL, quaCenterWeights, NULL, NULL, NULL, NULL}; 1853 PetscReal *sdWeights3[27] = {NULL, NULL, NULL, NULL, hexBottomWeights, NULL, NULL, NULL, NULL, NULL, hexFrontWeights, NULL, hexLeftWeights, hexCenterWeights, hexRightWeights, NULL, hexBackWeights, NULL, 1854 NULL, NULL, NULL, NULL, hexTopWeights, NULL, NULL, NULL, NULL}; 1855 PetscReal **sdWeights[4] = {NULL, NULL, sdWeights2, sdWeights3}; 1856 1857 /* Missing entries in serendipity cell, only works for 8-node quad and 20-node hex */ 1858 for (n = 0; n < elem->numNodes + s; ++n) { 1859 if (lexorder[n] >= 0) continue; 1860 for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] = 0.0; 1861 for (int bn = 0; bn < elem->numNodes + s; ++bn) { 1862 if (lexorder[bn] < 0) continue; 1863 const PetscReal *weights = sdWeights[coordDim][n]; 1864 const PetscInt bnode = elem->nodes[lexorder[bn]]; 1865 for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] += weights[bn] * (PetscReal)coords[bnode * 3 + d]; 1866 } 1867 } 1868 } 1869 PetscCall(DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES)); 1870 } 1871 PetscCall(PetscSectionDestroy(§ion)); 1872 PetscCall(PetscFree(cellCoords)); 1873 } else { 1874 PetscInt *nodeMap; 1875 double *coords = mesh ? mesh->nodelist->xyz : NULL; 1876 PetscScalar *pointCoords; 1877 1878 PetscCall(DMGetCoordinateSection(*dm, &cs)); 1879 PetscCall(PetscSectionSetNumFields(cs, 1)); 1880 PetscCall(PetscSectionSetFieldComponents(cs, 0, coordDim)); 1881 PetscCall(PetscSectionSetChart(cs, numCells, numCells + numVerts)); 1882 for (v = numCells; v < numCells + numVerts; ++v) { 1883 PetscCall(PetscSectionSetDof(cs, v, coordDim)); 1884 PetscCall(PetscSectionSetFieldDof(cs, v, 0, coordDim)); 1885 } 1886 PetscCall(PetscSectionSetUp(cs)); 1887 1888 /* We need to localize coordinates on cells */ 1889 if (periodic) { 1890 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)cdmCell), &csCell)); 1891 PetscCall(PetscSectionSetNumFields(csCell, 1)); 1892 PetscCall(PetscSectionSetFieldComponents(csCell, 0, coordDim)); 1893 PetscCall(PetscSectionSetChart(csCell, 0, numCells)); 1894 for (cell = 0; cell < numCells; ++cell) { 1895 if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) { 1896 GmshElement *elem = mesh->elements + cell; 1897 PetscInt dof = elem->numVerts * coordDim; 1898 1899 PetscCall(PetscSectionSetDof(csCell, cell, dof)); 1900 PetscCall(PetscSectionSetFieldDof(csCell, cell, 0, dof)); 1901 } 1902 } 1903 PetscCall(PetscSectionSetUp(csCell)); 1904 PetscCall(DMSetCellCoordinateSection(*dm, PETSC_DETERMINE, csCell)); 1905 } 1906 1907 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 1908 PetscCall(VecGetArray(coordinates, &pointCoords)); 1909 PetscCall(PetscMalloc1(numVerts, &nodeMap)); 1910 for (n = 0; n < numNodes; n++) 1911 if (mesh->vertexMap[n] >= 0) nodeMap[mesh->vertexMap[n]] = n; 1912 for (v = 0; v < numVerts; ++v) { 1913 PetscInt off, node = nodeMap[v]; 1914 1915 PetscCall(PetscSectionGetOffset(cs, numCells + v, &off)); 1916 for (d = 0; d < coordDim; ++d) pointCoords[off + d] = (PetscReal)coords[node * 3 + d]; 1917 } 1918 PetscCall(VecRestoreArray(coordinates, &pointCoords)); 1919 PetscCall(PetscFree(nodeMap)); 1920 1921 if (periodic) { 1922 PetscCall(DMCreateLocalVector(cdmCell, &coordinatesCell)); 1923 PetscCall(VecGetArray(coordinatesCell, &pointCoords)); 1924 for (cell = 0; cell < numCells; ++cell) { 1925 if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) { 1926 GmshElement *elem = mesh->elements + cell; 1927 PetscInt off, node; 1928 for (v = 0; v < elem->numVerts; ++v) cone[v] = elem->nodes[v]; 1929 PetscCall(DMPlexReorderCell(cdmCell, cell, cone)); 1930 PetscCall(PetscSectionGetOffset(csCell, cell, &off)); 1931 for (v = 0; v < elem->numVerts; ++v) 1932 for (node = cone[v], d = 0; d < coordDim; ++d) pointCoords[off++] = (PetscReal)coords[node * 3 + d]; 1933 } 1934 } 1935 PetscCall(VecSetBlockSize(coordinatesCell, coordDim)); 1936 PetscCall(VecRestoreArray(coordinatesCell, &pointCoords)); 1937 PetscCall(DMSetCellCoordinatesLocal(*dm, coordinatesCell)); 1938 PetscCall(VecDestroy(&coordinatesCell)); 1939 } 1940 PetscCall(PetscSectionDestroy(&csCell)); 1941 PetscCall(DMDestroy(&cdmCell)); 1942 } 1943 1944 PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates")); 1945 PetscCall(VecSetBlockSize(coordinates, coordDim)); 1946 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 1947 PetscCall(VecDestroy(&coordinates)); 1948 1949 PetscCall(GmshMeshDestroy(&mesh)); 1950 PetscCall(PetscBTDestroy(&periodicVerts)); 1951 PetscCall(PetscBTDestroy(&periodicCells)); 1952 1953 if (highOrder && project) { 1954 PetscFE fe; 1955 const char prefix[] = "dm_plex_gmsh_project_"; 1956 PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE; 1957 PetscDTNodeType nodeType = PETSCDTNODES_GAUSSJACOBI; 1958 1959 if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */ 1960 1961 PetscCall(GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe)); 1962 PetscCall(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view")); 1963 PetscCall(DMProjectCoordinates(*dm, fe)); 1964 PetscCall(PetscFEDestroy(&fe)); 1965 } 1966 1967 PetscCall(PetscLogEventEnd(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL)); 1968 PetscFunctionReturn(0); 1969 } 1970