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