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 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 PetscCheck(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 PetscCheck(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 PetscCheck(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 PetscCheck(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 PetscCheck(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 PetscCheck(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 PetscCheck(snum == 3,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line); 1066 PetscCheck(version >= 2.2,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version); 1067 PetscCheck((int)version != 3,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version); 1068 PetscCheck(version <= 4.1,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version); 1069 PetscCheck(!gmsh->binary || fileType,PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is binary but Gmsh file is ASCII"); 1070 PetscCheck(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 PetscCheck(fileFormat > 40 || dataSize == sizeof(double),PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Data size %d is not valid for a Gmsh file", dataSize); 1073 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); 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 PetscCheck(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 PetscCheck(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 PetscCheck(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 PetscCheck(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 PetscCheck(gmsh->nodeMap[tag] < 0,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Repeated node tag %" PetscInt_FMT, 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%" PetscInt_FMT, 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 PetscCheck(snum == 2,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line); 1405 PetscCheck(version >= 2.2,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version); 1406 PetscCheck((int)version != 3,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version); 1407 PetscCheck(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 1461 PetscFunctionBegin; 1462 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1463 PetscObjectOptionsBegin((PetscObject)viewer); 1464 PetscOptionsHeadBegin(PetscOptionsObject,"DMPlex Gmsh options"); 1465 PetscCall(PetscOptionsBool("-dm_plex_gmsh_hybrid", "Generate hybrid cell bounds", "DMPlexCreateGmsh", hybrid, &hybrid, NULL)); 1466 PetscCall(PetscOptionsBool("-dm_plex_gmsh_periodic","Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL)); 1467 PetscCall(PetscOptionsBool("-dm_plex_gmsh_highorder","Generate high-order coordinates", "DMPlexCreateGmsh", highOrder, &highOrder, &highOrderSet)); 1468 PetscCall(PetscOptionsBool("-dm_plex_gmsh_project", "Project high-order coordinates to a different space", "DMPlexCreateGmsh", project, &project, NULL)); 1469 PetscCall(PetscOptionsBool("-dm_plex_gmsh_use_marker", "Generate marker label", "DMPlexCreateGmsh", usemarker, &usemarker, NULL)); 1470 PetscCall(PetscOptionsBool("-dm_plex_gmsh_use_regions", "Generate labels with region names", "DMPlexCreateGmsh", useregions, &useregions, NULL)); 1471 PetscCall(PetscOptionsBool("-dm_plex_gmsh_mark_vertices", "Add vertices to generated labels", "DMPlexCreateGmsh", markvertices, &markvertices, NULL)); 1472 PetscCall(PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE)); 1473 PetscOptionsHeadEnd(); 1474 PetscOptionsEnd(); 1475 1476 PetscCall(GmshCellInfoSetUp()); 1477 1478 PetscCall(DMCreate(comm, dm)); 1479 PetscCall(DMSetType(*dm, DMPLEX)); 1480 PetscCall(PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,NULL,NULL,NULL)); 1481 1482 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary)); 1483 1484 /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */ 1485 if (binary) { 1486 parentviewer = viewer; 1487 PetscCall(PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer)); 1488 } 1489 1490 if (rank == 0) { 1491 GmshFile gmsh[1]; 1492 char line[PETSC_MAX_PATH_LEN]; 1493 PetscBool match; 1494 1495 PetscCall(PetscArrayzero(gmsh,1)); 1496 gmsh->viewer = viewer; 1497 gmsh->binary = binary; 1498 1499 PetscCall(GmshMeshCreate(&mesh)); 1500 1501 /* Read mesh format */ 1502 PetscCall(GmshReadSection(gmsh, line)); 1503 PetscCall(GmshExpect(gmsh, "$MeshFormat", line)); 1504 PetscCall(GmshReadMeshFormat(gmsh)); 1505 PetscCall(GmshReadEndSection(gmsh, "$EndMeshFormat", line)); 1506 1507 /* OPTIONAL Read physical names */ 1508 PetscCall(GmshReadSection(gmsh, line)); 1509 PetscCall(GmshMatch(gmsh, "$PhysicalNames", line, &match)); 1510 if (match) { 1511 PetscCall(GmshExpect(gmsh, "$PhysicalNames", line)); 1512 PetscCall(GmshReadPhysicalNames(gmsh, mesh)); 1513 PetscCall(GmshReadEndSection(gmsh, "$EndPhysicalNames", line)); 1514 /* Initial read for entity section */ 1515 PetscCall(GmshReadSection(gmsh, line)); 1516 } 1517 1518 /* Read entities */ 1519 if (gmsh->fileFormat >= 40) { 1520 PetscCall(GmshExpect(gmsh, "$Entities", line)); 1521 PetscCall(GmshReadEntities(gmsh, mesh)); 1522 PetscCall(GmshReadEndSection(gmsh, "$EndEntities", line)); 1523 /* Initial read for nodes section */ 1524 PetscCall(GmshReadSection(gmsh, line)); 1525 } 1526 1527 /* Read nodes */ 1528 PetscCall(GmshExpect(gmsh, "$Nodes", line)); 1529 PetscCall(GmshReadNodes(gmsh, mesh)); 1530 PetscCall(GmshReadEndSection(gmsh, "$EndNodes", line)); 1531 1532 /* Read elements */ 1533 PetscCall(GmshReadSection(gmsh, line)); 1534 PetscCall(GmshExpect(gmsh, "$Elements", line)); 1535 PetscCall(GmshReadElements(gmsh, mesh)); 1536 PetscCall(GmshReadEndSection(gmsh, "$EndElements", line)); 1537 1538 /* Read periodic section (OPTIONAL) */ 1539 if (periodic) { 1540 PetscCall(GmshReadSection(gmsh, line)); 1541 PetscCall(GmshMatch(gmsh, "$Periodic", line, &periodic)); 1542 } 1543 if (periodic) { 1544 PetscCall(GmshExpect(gmsh, "$Periodic", line)); 1545 PetscCall(GmshReadPeriodic(gmsh, mesh)); 1546 PetscCall(GmshReadEndSection(gmsh, "$EndPeriodic", line)); 1547 } 1548 1549 PetscCall(PetscFree(gmsh->wbuf)); 1550 PetscCall(PetscFree(gmsh->sbuf)); 1551 PetscCall(PetscFree(gmsh->nbuf)); 1552 1553 dim = mesh->dim; 1554 order = mesh->order; 1555 numNodes = mesh->numNodes; 1556 numElems = mesh->numElems; 1557 numVerts = mesh->numVerts; 1558 numCells = mesh->numCells; 1559 1560 { 1561 GmshElement *elemA = mesh->numCells > 0 ? mesh->elements : NULL; 1562 GmshElement *elemB = elemA ? elemA + mesh->numCells - 1 : NULL; 1563 int ptA = elemA ? GmshCellMap[elemA->cellType].polytope : -1; 1564 int ptB = elemB ? GmshCellMap[elemB->cellType].polytope : -1; 1565 isSimplex = (ptA == GMSH_QUA || ptA == GMSH_HEX) ? PETSC_FALSE : PETSC_TRUE; 1566 isHybrid = (ptA == ptB) ? PETSC_FALSE : PETSC_TRUE; 1567 hasTetra = (ptA == GMSH_TET) ? PETSC_TRUE : PETSC_FALSE; 1568 } 1569 } 1570 1571 if (parentviewer) { 1572 PetscCall(PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer)); 1573 } 1574 1575 { 1576 int buf[6]; 1577 1578 buf[0] = (int)dim; 1579 buf[1] = (int)order; 1580 buf[2] = periodic; 1581 buf[3] = isSimplex; 1582 buf[4] = isHybrid; 1583 buf[5] = hasTetra; 1584 1585 PetscCallMPI(MPI_Bcast(buf, 6, MPI_INT, 0, comm)); 1586 1587 dim = buf[0]; 1588 order = buf[1]; 1589 periodic = buf[2] ? PETSC_TRUE : PETSC_FALSE; 1590 isSimplex = buf[3] ? PETSC_TRUE : PETSC_FALSE; 1591 isHybrid = buf[4] ? PETSC_TRUE : PETSC_FALSE; 1592 hasTetra = buf[5] ? PETSC_TRUE : PETSC_FALSE; 1593 } 1594 1595 if (!highOrderSet) highOrder = (order > 1) ? PETSC_TRUE : PETSC_FALSE; 1596 PetscCheck(!highOrder || !isHybrid,comm, PETSC_ERR_SUP, "No support for discretization on hybrid meshes yet"); 1597 1598 /* We do not want this label automatically computed, instead we fill it here */ 1599 PetscCall(DMCreateLabel(*dm, "celltype")); 1600 1601 /* Allocate the cell-vertex mesh */ 1602 PetscCall(DMPlexSetChart(*dm, 0, numCells+numVerts)); 1603 for (cell = 0; cell < numCells; ++cell) { 1604 GmshElement *elem = mesh->elements + cell; 1605 DMPolytopeType ctype = DMPolytopeTypeFromGmsh(elem->cellType); 1606 if (hybrid && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR; 1607 PetscCall(DMPlexSetConeSize(*dm, cell, elem->numVerts)); 1608 PetscCall(DMPlexSetCellType(*dm, cell, ctype)); 1609 } 1610 for (v = numCells; v < numCells+numVerts; ++v) { 1611 PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT)); 1612 } 1613 PetscCall(DMSetUp(*dm)); 1614 1615 /* Add cell-vertex connections */ 1616 for (cell = 0; cell < numCells; ++cell) { 1617 GmshElement *elem = mesh->elements + cell; 1618 for (v = 0; v < elem->numVerts; ++v) { 1619 const PetscInt nn = elem->nodes[v]; 1620 const PetscInt vv = mesh->vertexMap[nn]; 1621 cone[v] = numCells + vv; 1622 } 1623 PetscCall(DMPlexReorderCell(*dm, cell, cone)); 1624 PetscCall(DMPlexSetCone(*dm, cell, cone)); 1625 } 1626 1627 PetscCall(DMSetDimension(*dm, dim)); 1628 PetscCall(DMPlexSymmetrize(*dm)); 1629 PetscCall(DMPlexStratify(*dm)); 1630 if (interpolate) { 1631 DM idm; 1632 1633 PetscCall(DMPlexInterpolate(*dm, &idm)); 1634 PetscCall(DMDestroy(dm)); 1635 *dm = idm; 1636 } 1637 1638 /* Create the label "marker" over the whole boundary */ 1639 PetscCheck(!usemarker || interpolate || dim <= 1,comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation"); 1640 if (rank == 0 && usemarker) { 1641 PetscInt f, fStart, fEnd; 1642 1643 PetscCall(DMCreateLabel(*dm, "marker")); 1644 PetscCall(DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd)); 1645 for (f = fStart; f < fEnd; ++f) { 1646 PetscInt suppSize; 1647 1648 PetscCall(DMPlexGetSupportSize(*dm, f, &suppSize)); 1649 if (suppSize == 1) { 1650 PetscInt *cone = NULL, coneSize, p; 1651 1652 PetscCall(DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone)); 1653 for (p = 0; p < coneSize; p += 2) { 1654 PetscCall(DMSetLabelValue_Fast(*dm, &marker, "marker", cone[p], 1)); 1655 } 1656 PetscCall(DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone)); 1657 } 1658 } 1659 } 1660 1661 if (rank == 0) { 1662 const PetscInt Nr = useregions ? mesh->numRegions : 0; 1663 PetscInt vStart, vEnd; 1664 1665 PetscCall(PetscCalloc1(Nr, ®ionSets)); 1666 PetscCall(DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd)); 1667 for (cell = 0, e = 0; e < numElems; ++e) { 1668 GmshElement *elem = mesh->elements + e; 1669 1670 /* Create cell sets */ 1671 if (elem->dim == dim && dim > 0) { 1672 if (elem->numTags > 0) { 1673 const PetscInt tag = elem->tags[0]; 1674 PetscInt r; 1675 1676 if (!Nr) PetscCall(DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", cell, tag)); 1677 for (r = 0; r < Nr; ++r) { 1678 if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], cell, tag)); 1679 } 1680 } 1681 cell++; 1682 } 1683 1684 /* Create face sets */ 1685 if (interpolate && elem->dim == dim-1) { 1686 PetscInt joinSize; 1687 const PetscInt *join = NULL; 1688 const PetscInt tag = elem->tags[0]; 1689 PetscInt r; 1690 1691 /* Find the relevant facet with vertex joins */ 1692 for (v = 0; v < elem->numVerts; ++v) { 1693 const PetscInt nn = elem->nodes[v]; 1694 const PetscInt vv = mesh->vertexMap[nn]; 1695 cone[v] = vStart + vv; 1696 } 1697 PetscCall(DMPlexGetFullJoin(*dm, elem->numVerts, cone, &joinSize, &join)); 1698 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); 1699 if (!Nr) PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], tag)); 1700 for (r = 0; r < Nr; ++r) { 1701 if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], join[0], tag)); 1702 } 1703 PetscCall(DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join)); 1704 } 1705 1706 /* Create vertex sets */ 1707 if (elem->dim == 0) { 1708 if (elem->numTags > 0) { 1709 const PetscInt nn = elem->nodes[0]; 1710 const PetscInt vv = mesh->vertexMap[nn]; 1711 const PetscInt tag = elem->tags[0]; 1712 PetscInt r; 1713 1714 if (!Nr) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag)); 1715 for (r = 0; r < Nr; ++r) { 1716 if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag)); 1717 } 1718 } 1719 } 1720 } 1721 if (markvertices) { 1722 for (v = 0; v < numNodes; ++v) { 1723 const PetscInt vv = mesh->vertexMap[v]; 1724 const PetscInt tag = mesh->nodelist->tag[v]; 1725 PetscInt r; 1726 1727 if (tag != -1) { 1728 if (!Nr) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag)); 1729 for (r = 0; r < Nr; ++r) { 1730 if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag)); 1731 } 1732 } 1733 } 1734 } 1735 PetscCall(PetscFree(regionSets)); 1736 } 1737 1738 { /* Create Cell/Face/Vertex Sets labels at all processes */ 1739 enum {n = 4}; 1740 PetscBool flag[n]; 1741 1742 flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE; 1743 flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE; 1744 flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE; 1745 flag[3] = marker ? PETSC_TRUE : PETSC_FALSE; 1746 PetscCallMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm)); 1747 if (flag[0]) PetscCall(DMCreateLabel(*dm, "Cell Sets")); 1748 if (flag[1]) PetscCall(DMCreateLabel(*dm, "Face Sets")); 1749 if (flag[2]) PetscCall(DMCreateLabel(*dm, "Vertex Sets")); 1750 if (flag[3]) PetscCall(DMCreateLabel(*dm, "marker")); 1751 } 1752 1753 if (periodic) { 1754 PetscCall(PetscBTCreate(numVerts, &periodicVerts)); 1755 for (n = 0; n < numNodes; ++n) { 1756 if (mesh->vertexMap[n] >= 0) { 1757 if (PetscUnlikely(mesh->periodMap[n] != n)) { 1758 PetscInt m = mesh->periodMap[n]; 1759 PetscCall(PetscBTSet(periodicVerts, mesh->vertexMap[n])); 1760 PetscCall(PetscBTSet(periodicVerts, mesh->vertexMap[m])); 1761 } 1762 } 1763 } 1764 PetscCall(PetscBTCreate(numCells, &periodicCells)); 1765 for (cell = 0; cell < numCells; ++cell) { 1766 GmshElement *elem = mesh->elements + cell; 1767 for (v = 0; v < elem->numVerts; ++v) { 1768 PetscInt nn = elem->nodes[v]; 1769 PetscInt vv = mesh->vertexMap[nn]; 1770 if (PetscUnlikely(PetscBTLookup(periodicVerts, vv))) { 1771 PetscCall(PetscBTSet(periodicCells, cell)); break; 1772 } 1773 } 1774 } 1775 } 1776 1777 /* Setup coordinate DM */ 1778 if (coordDim < 0) coordDim = dim; 1779 PetscCall(DMSetCoordinateDim(*dm, coordDim)); 1780 PetscCall(DMGetCoordinateDM(*dm, &cdm)); 1781 if (highOrder) { 1782 PetscFE fe; 1783 PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE; 1784 PetscDTNodeType nodeType = PETSCDTNODES_EQUISPACED; 1785 1786 if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */ 1787 1788 PetscCall(GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe)); 1789 PetscCall(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view")); 1790 PetscCall(DMSetField(cdm, 0, NULL, (PetscObject) fe)); 1791 PetscCall(PetscFEDestroy(&fe)); 1792 PetscCall(DMCreateDS(cdm)); 1793 } 1794 1795 /* Create coordinates */ 1796 if (highOrder) { 1797 1798 PetscInt maxDof = GmshNumNodes_HEX(order)*coordDim; 1799 double *coords = mesh ? mesh->nodelist->xyz : NULL; 1800 PetscSection section; 1801 PetscScalar *cellCoords; 1802 1803 PetscCall(DMSetLocalSection(cdm, NULL)); 1804 PetscCall(DMGetLocalSection(cdm, &coordSection)); 1805 PetscCall(PetscSectionClone(coordSection, §ion)); 1806 PetscCall(DMPlexSetClosurePermutationTensor(cdm, 0, section)); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */ 1807 1808 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 1809 PetscCall(PetscMalloc1(maxDof, &cellCoords)); 1810 for (cell = 0; cell < numCells; ++cell) { 1811 GmshElement *elem = mesh->elements + cell; 1812 const int *lexorder = GmshCellMap[elem->cellType].lexorder(); 1813 int s = 0; 1814 for (n = 0; n < elem->numNodes; ++n) { 1815 while (lexorder[n+s] < 0) ++s; 1816 const PetscInt node = elem->nodes[lexorder[n+s]]; 1817 for (d = 0; d < coordDim; ++d) cellCoords[(n+s)*coordDim+d] = (PetscReal) coords[node*3+d]; 1818 } 1819 if (s) { 1820 /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */ 1821 PetscReal quaCenterWeights[9] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25}; 1822 /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */ 1823 PetscReal hexBottomWeights[27] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25, 1824 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1825 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; 1826 PetscReal hexFrontWeights[27] = {-0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1827 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1828 -0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; 1829 PetscReal hexLeftWeights[27] = {-0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 1830 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 1831 -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0}; 1832 PetscReal hexRightWeights[27] = { 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 1833 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1834 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25}; 1835 PetscReal hexBackWeights[27] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25, 1836 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 1837 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25}; 1838 PetscReal hexTopWeights[27] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1839 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1840 -0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25}; 1841 PetscReal hexCenterWeights[27] = {-0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25, 1842 0.25, 0.0, 0.25, 0.0, 0.0, 0.0, 0.25, 0.0, 0.25, 1843 -0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25}; 1844 PetscReal *sdWeights2[9] = {NULL, NULL, NULL, NULL, quaCenterWeights, NULL, NULL, NULL, NULL}; 1845 PetscReal *sdWeights3[27] = {NULL, NULL, NULL, NULL, hexBottomWeights, NULL, NULL, NULL, NULL, 1846 NULL, hexFrontWeights, NULL, hexLeftWeights, hexCenterWeights, hexRightWeights, NULL, hexBackWeights, NULL, 1847 NULL, NULL, NULL, NULL, hexTopWeights, NULL, NULL, NULL, NULL}; 1848 PetscReal **sdWeights[4] = {NULL, NULL, sdWeights2, sdWeights3}; 1849 1850 /* Missing entries in serendipity cell, only works for 8-node quad and 20-node hex */ 1851 for (n = 0; n < elem->numNodes+s; ++n) { 1852 if (lexorder[n] >= 0) continue; 1853 for (d = 0; d < coordDim; ++d) cellCoords[n*coordDim+d] = 0.0; 1854 for (int bn = 0; bn < elem->numNodes+s; ++bn) { 1855 if (lexorder[bn] < 0) continue; 1856 const PetscReal *weights = sdWeights[coordDim][n]; 1857 const PetscInt bnode = elem->nodes[lexorder[bn]]; 1858 for (d = 0; d < coordDim; ++d) cellCoords[n*coordDim+d] += weights[bn] * (PetscReal) coords[bnode*3+d]; 1859 } 1860 } 1861 } 1862 PetscCall(DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES)); 1863 } 1864 PetscCall(PetscSectionDestroy(§ion)); 1865 PetscCall(PetscFree(cellCoords)); 1866 1867 } else { 1868 1869 PetscInt *nodeMap; 1870 double *coords = mesh ? mesh->nodelist->xyz : NULL; 1871 PetscScalar *pointCoords; 1872 1873 PetscCall(DMGetLocalSection(cdm, &coordSection)); 1874 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 1875 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, coordDim)); 1876 if (periodic) { /* we need to localize coordinates on cells */ 1877 PetscCall(PetscSectionSetChart(coordSection, 0, numCells+numVerts)); 1878 } else { 1879 PetscCall(PetscSectionSetChart(coordSection, numCells, numCells+numVerts)); 1880 } 1881 if (periodic) { 1882 for (cell = 0; cell < numCells; ++cell) { 1883 if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) { 1884 GmshElement *elem = mesh->elements + cell; 1885 PetscInt dof = elem->numVerts * coordDim; 1886 PetscCall(PetscSectionSetDof(coordSection, cell, dof)); 1887 PetscCall(PetscSectionSetFieldDof(coordSection, cell, 0, dof)); 1888 } 1889 } 1890 } 1891 for (v = numCells; v < numCells+numVerts; ++v) { 1892 PetscCall(PetscSectionSetDof(coordSection, v, coordDim)); 1893 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, coordDim)); 1894 } 1895 PetscCall(PetscSectionSetUp(coordSection)); 1896 1897 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 1898 PetscCall(VecGetArray(coordinates, &pointCoords)); 1899 if (periodic) { 1900 for (cell = 0; cell < numCells; ++cell) { 1901 if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) { 1902 GmshElement *elem = mesh->elements + cell; 1903 PetscInt off, node; 1904 for (v = 0; v < elem->numVerts; ++v) 1905 cone[v] = elem->nodes[v]; 1906 PetscCall(DMPlexReorderCell(cdm, cell, cone)); 1907 PetscCall(PetscSectionGetOffset(coordSection, cell, &off)); 1908 for (v = 0; v < elem->numVerts; ++v) 1909 for (node = cone[v], d = 0; d < coordDim; ++d) 1910 pointCoords[off++] = (PetscReal) coords[node*3+d]; 1911 } 1912 } 1913 } 1914 PetscCall(PetscMalloc1(numVerts, &nodeMap)); 1915 for (n = 0; n < numNodes; n++) 1916 if (mesh->vertexMap[n] >= 0) 1917 nodeMap[mesh->vertexMap[n]] = n; 1918 for (v = 0; v < numVerts; ++v) { 1919 PetscInt off, node = nodeMap[v]; 1920 PetscCall(PetscSectionGetOffset(coordSection, numCells + v, &off)); 1921 for (d = 0; d < coordDim; ++d) 1922 pointCoords[off+d] = (PetscReal) coords[node*3+d]; 1923 } 1924 PetscCall(PetscFree(nodeMap)); 1925 PetscCall(VecRestoreArray(coordinates, &pointCoords)); 1926 1927 } 1928 1929 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 1930 PetscCall(VecSetBlockSize(coordinates, coordDim)); 1931 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 1932 PetscCall(VecDestroy(&coordinates)); 1933 PetscCall(DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL)); 1934 1935 PetscCall(GmshMeshDestroy(&mesh)); 1936 PetscCall(PetscBTDestroy(&periodicVerts)); 1937 PetscCall(PetscBTDestroy(&periodicCells)); 1938 1939 if (highOrder && project) { 1940 PetscFE fe; 1941 const char prefix[] = "dm_plex_gmsh_project_"; 1942 PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE; 1943 PetscDTNodeType nodeType = PETSCDTNODES_GAUSSJACOBI; 1944 1945 if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */ 1946 1947 PetscCall(GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe)); 1948 PetscCall(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view")); 1949 PetscCall(DMProjectCoordinates(*dm, fe)); 1950 PetscCall(PetscFEDestroy(&fe)); 1951 } 1952 1953 PetscCall(PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,NULL,NULL,NULL)); 1954 PetscFunctionReturn(0); 1955 } 1956