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