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