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