1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 #include <petsc/private/hashmapi.h> 4 5 #include <../src/dm/impls/plex/gmshlex.h> 6 7 #define GMSH_LEXORDER_ITEM(T, p) \ 8 static int *GmshLexOrder_##T##_##p(void) \ 9 { \ 10 static int Gmsh_LexOrder_##T##_##p[GmshNumNodes_##T(p)] = {-1}; \ 11 int *lex = Gmsh_LexOrder_##T##_##p; \ 12 if (lex[0] == -1) (void)GmshLexOrder_##T(p, lex, 0); \ 13 return lex; \ 14 } 15 16 static int *GmshLexOrder_QUA_2_Serendipity(void) 17 { 18 static int Gmsh_LexOrder_QUA_2_Serendipity[9] = {-1}; 19 int *lex = Gmsh_LexOrder_QUA_2_Serendipity; 20 if (lex[0] == -1) { 21 /* Vertices */ 22 lex[0] = 0; lex[2] = 1; lex[8] = 2; lex[6] = 3; 23 /* Edges */ 24 lex[1] = 4; lex[5] = 5; lex[7] = 6; lex[3] = 7; 25 /* Cell */ 26 lex[4] = -8; 27 } 28 return lex; 29 } 30 31 static int *GmshLexOrder_HEX_2_Serendipity(void) 32 { 33 static int Gmsh_LexOrder_HEX_2_Serendipity[27] = {-1}; 34 int *lex = Gmsh_LexOrder_HEX_2_Serendipity; 35 if (lex[0] == -1) { 36 /* Vertices */ 37 lex[ 0] = 0; lex[ 2] = 1; lex[ 8] = 2; lex[ 6] = 3; 38 lex[18] = 4; lex[20] = 5; lex[26] = 6; lex[24] = 7; 39 /* Edges */ 40 lex[ 1] = 8; lex[ 3] = 9; lex[ 9] = 10; lex[ 5] = 11; 41 lex[11] = 12; lex[ 7] = 13; lex[17] = 14; lex[15] = 15; 42 lex[19] = 16; lex[21] = 17; lex[23] = 18; lex[25] = 19; 43 /* Faces */ 44 lex[ 4] = -20; lex[10] = -21; lex[12] = -22; lex[14] = -23; 45 lex[16] = -24; lex[22] = -25; 46 /* Cell */ 47 lex[13] = -26; 48 } 49 return lex; 50 } 51 52 #define GMSH_LEXORDER_LIST(T) \ 53 GMSH_LEXORDER_ITEM(T, 1) \ 54 GMSH_LEXORDER_ITEM(T, 2) \ 55 GMSH_LEXORDER_ITEM(T, 3) \ 56 GMSH_LEXORDER_ITEM(T, 4) \ 57 GMSH_LEXORDER_ITEM(T, 5) \ 58 GMSH_LEXORDER_ITEM(T, 6) \ 59 GMSH_LEXORDER_ITEM(T, 7) \ 60 GMSH_LEXORDER_ITEM(T, 8) \ 61 GMSH_LEXORDER_ITEM(T, 9) \ 62 GMSH_LEXORDER_ITEM(T, 10) 63 64 GMSH_LEXORDER_ITEM(VTX, 0) 65 GMSH_LEXORDER_LIST(SEG) 66 GMSH_LEXORDER_LIST(TRI) 67 GMSH_LEXORDER_LIST(QUA) 68 GMSH_LEXORDER_LIST(TET) 69 GMSH_LEXORDER_LIST(HEX) 70 GMSH_LEXORDER_LIST(PRI) 71 GMSH_LEXORDER_LIST(PYR) 72 73 typedef enum { 74 GMSH_VTX = 0, 75 GMSH_SEG = 1, 76 GMSH_TRI = 2, 77 GMSH_QUA = 3, 78 GMSH_TET = 4, 79 GMSH_HEX = 5, 80 GMSH_PRI = 6, 81 GMSH_PYR = 7, 82 GMSH_NUM_POLYTOPES = 8 83 } GmshPolytopeType; 84 85 typedef struct { 86 int cellType; 87 int polytope; 88 int dim; 89 int order; 90 int numVerts; 91 int numNodes; 92 int* (*lexorder)(void); 93 } GmshCellInfo; 94 95 #define GmshCellEntry(cellType, polytope, dim, order) \ 96 {cellType, GMSH_##polytope, dim, order, \ 97 GmshNumNodes_##polytope(1), \ 98 GmshNumNodes_##polytope(order), \ 99 GmshLexOrder_##polytope##_##order} 100 101 static const GmshCellInfo GmshCellTable[] = { 102 GmshCellEntry( 15, VTX, 0, 0), 103 104 GmshCellEntry( 1, SEG, 1, 1), 105 GmshCellEntry( 8, SEG, 1, 2), 106 GmshCellEntry( 26, SEG, 1, 3), 107 GmshCellEntry( 27, SEG, 1, 4), 108 GmshCellEntry( 28, SEG, 1, 5), 109 GmshCellEntry( 62, SEG, 1, 6), 110 GmshCellEntry( 63, SEG, 1, 7), 111 GmshCellEntry( 64, SEG, 1, 8), 112 GmshCellEntry( 65, SEG, 1, 9), 113 GmshCellEntry( 66, SEG, 1, 10), 114 115 GmshCellEntry( 2, TRI, 2, 1), 116 GmshCellEntry( 9, TRI, 2, 2), 117 GmshCellEntry( 21, TRI, 2, 3), 118 GmshCellEntry( 23, TRI, 2, 4), 119 GmshCellEntry( 25, TRI, 2, 5), 120 GmshCellEntry( 42, TRI, 2, 6), 121 GmshCellEntry( 43, TRI, 2, 7), 122 GmshCellEntry( 44, TRI, 2, 8), 123 GmshCellEntry( 45, TRI, 2, 9), 124 GmshCellEntry( 46, TRI, 2, 10), 125 126 GmshCellEntry( 3, QUA, 2, 1), 127 GmshCellEntry( 10, QUA, 2, 2), 128 {16, GMSH_QUA, 2, 2, 4, 8, GmshLexOrder_QUA_2_Serendipity}, 129 GmshCellEntry( 36, QUA, 2, 3), 130 GmshCellEntry( 37, QUA, 2, 4), 131 GmshCellEntry( 38, QUA, 2, 5), 132 GmshCellEntry( 47, QUA, 2, 6), 133 GmshCellEntry( 48, QUA, 2, 7), 134 GmshCellEntry( 49, QUA, 2, 8), 135 GmshCellEntry( 50, QUA, 2, 9), 136 GmshCellEntry( 51, QUA, 2, 10), 137 138 GmshCellEntry( 4, TET, 3, 1), 139 GmshCellEntry( 11, TET, 3, 2), 140 GmshCellEntry( 29, TET, 3, 3), 141 GmshCellEntry( 30, TET, 3, 4), 142 GmshCellEntry( 31, TET, 3, 5), 143 GmshCellEntry( 71, TET, 3, 6), 144 GmshCellEntry( 72, TET, 3, 7), 145 GmshCellEntry( 73, TET, 3, 8), 146 GmshCellEntry( 74, TET, 3, 9), 147 GmshCellEntry( 75, TET, 3, 10), 148 149 GmshCellEntry( 5, HEX, 3, 1), 150 GmshCellEntry( 12, HEX, 3, 2), 151 {17, GMSH_HEX, 3, 2, 8, 20, GmshLexOrder_HEX_2_Serendipity}, 152 GmshCellEntry( 92, HEX, 3, 3), 153 GmshCellEntry( 93, HEX, 3, 4), 154 GmshCellEntry( 94, HEX, 3, 5), 155 GmshCellEntry( 95, HEX, 3, 6), 156 GmshCellEntry( 96, HEX, 3, 7), 157 GmshCellEntry( 97, HEX, 3, 8), 158 GmshCellEntry( 98, HEX, 3, 9), 159 GmshCellEntry( -1, HEX, 3, 10), 160 161 GmshCellEntry( 6, PRI, 3, 1), 162 GmshCellEntry( 13, PRI, 3, 2), 163 GmshCellEntry( 90, PRI, 3, 3), 164 GmshCellEntry( 91, PRI, 3, 4), 165 GmshCellEntry(106, PRI, 3, 5), 166 GmshCellEntry(107, PRI, 3, 6), 167 GmshCellEntry(108, PRI, 3, 7), 168 GmshCellEntry(109, PRI, 3, 8), 169 GmshCellEntry(110, PRI, 3, 9), 170 GmshCellEntry( -1, PRI, 3, 10), 171 172 GmshCellEntry( 7, PYR, 3, 1), 173 GmshCellEntry( 14, PYR, 3, 2), 174 GmshCellEntry(118, PYR, 3, 3), 175 GmshCellEntry(119, PYR, 3, 4), 176 GmshCellEntry(120, PYR, 3, 5), 177 GmshCellEntry(121, PYR, 3, 6), 178 GmshCellEntry(122, PYR, 3, 7), 179 GmshCellEntry(123, PYR, 3, 8), 180 GmshCellEntry(124, PYR, 3, 9), 181 GmshCellEntry( -1, PYR, 3, 10) 182 183 #if 0 184 {20, GMSH_TRI, 2, 3, 3, 9, NULL}, 185 {18, GMSH_PRI, 3, 2, 6, 15, NULL}, 186 {19, GMSH_PYR, 3, 2, 5, 13, NULL}, 187 #endif 188 }; 189 190 static GmshCellInfo GmshCellMap[150]; 191 192 static PetscErrorCode GmshCellInfoSetUp(void) 193 { 194 size_t i, n; 195 static PetscBool called = PETSC_FALSE; 196 197 if (called) return 0; 198 PetscFunctionBegin; 199 called = PETSC_TRUE; 200 n = sizeof(GmshCellMap)/sizeof(GmshCellMap[0]); 201 for (i = 0; i < n; ++i) { 202 GmshCellMap[i].cellType = -1; 203 GmshCellMap[i].polytope = -1; 204 } 205 n = sizeof(GmshCellTable)/sizeof(GmshCellTable[0]); 206 for (i = 0; i < n; ++i) { 207 if (GmshCellTable[i].cellType <= 0) continue; 208 GmshCellMap[GmshCellTable[i].cellType] = GmshCellTable[i]; 209 } 210 PetscFunctionReturn(0); 211 } 212 213 #define GmshCellTypeCheck(ct) PetscMacroReturnStandard( \ 214 const int _ct_ = (int)ct; \ 215 PetscCheck(_ct_ >= 0 && _ct_ < (int)(sizeof(GmshCellMap)/sizeof(GmshCellMap[0])), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid Gmsh element type %d", _ct_); \ 216 PetscCheck(GmshCellMap[_ct_].cellType == _ct_, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported Gmsh element type %d", _ct_); \ 217 PetscCheck(GmshCellMap[_ct_].polytope != -1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported Gmsh element type %d", _ct_); \ 218 ) 219 220 typedef struct { 221 PetscViewer viewer; 222 int fileFormat; 223 int dataSize; 224 PetscBool binary; 225 PetscBool byteSwap; 226 size_t wlen; 227 void *wbuf; 228 size_t slen; 229 void *sbuf; 230 PetscInt *nbuf; 231 PetscInt nodeStart; 232 PetscInt nodeEnd; 233 PetscInt *nodeMap; 234 } GmshFile; 235 236 static PetscErrorCode GmshBufferGet(GmshFile *gmsh, size_t count, size_t eltsize, void *buf) 237 { 238 size_t size = count * eltsize; 239 240 PetscFunctionBegin; 241 if (gmsh->wlen < size) { 242 CHKERRQ(PetscFree(gmsh->wbuf)); 243 CHKERRQ(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 CHKERRQ(PetscFree(gmsh->sbuf)); 258 CHKERRQ(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 CHKERRQ(PetscViewerRead(gmsh->viewer, buf, count, NULL, dtype)); 269 if (gmsh->byteSwap) CHKERRQ(PetscByteSwap(buf, dtype, count)); 270 PetscFunctionReturn(0); 271 } 272 273 static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscInt count) 274 { 275 PetscFunctionBegin; 276 CHKERRQ(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 CHKERRQ(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 CHKERRQ(GmshMatch(gmsh, Section, line, &match)); 293 PetscCheckFalse(!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 CHKERRQ(GmshReadString(gmsh, line, 1)); 304 CHKERRQ(GmshMatch(gmsh, "$Comments", line, &match)); 305 if (!match) break; 306 while (PETSC_TRUE) { 307 CHKERRQ(GmshReadString(gmsh, line, 1)); 308 CHKERRQ(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 CHKERRQ(GmshReadString(gmsh, line, 1)); 319 CHKERRQ(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 CHKERRQ(GmshRead(gmsh, buf, count, PETSC_INT)); 331 } else if (dataSize == sizeof(int)) { 332 int *ibuf = NULL; 333 CHKERRQ(GmshBufferSizeGet(gmsh, count, &ibuf)); 334 CHKERRQ(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 CHKERRQ(GmshBufferSizeGet(gmsh, count, &ibuf)); 339 CHKERRQ(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 CHKERRQ(GmshBufferSizeGet(gmsh, count, &ibuf)); 344 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscNew(entities)); 383 for (dim = 0; dim < 4; ++dim) { 384 CHKERRQ(PetscCalloc1(count[dim], &(*entities)->entity[dim])); 385 CHKERRQ(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 CHKERRQ(PetscFree((*entities)->entity[dim])); 398 CHKERRQ(PetscHMapIDestroy(&(*entities)->entityMap[dim])); 399 } 400 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscNew(nodes)); 434 CHKERRQ(PetscMalloc1(count*1, &(*nodes)->id)); 435 CHKERRQ(PetscMalloc1(count*3, &(*nodes)->xyz)); 436 CHKERRQ(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 CHKERRQ(PetscFree((*nodes)->id)); 445 CHKERRQ(PetscFree((*nodes)->xyz)); 446 CHKERRQ(PetscFree((*nodes)->tag)); 447 CHKERRQ(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 CHKERRQ(PetscCalloc1(count, elements)); 466 PetscFunctionReturn(0); 467 } 468 469 static PetscErrorCode GmshElementsDestroy(GmshElement **elements) 470 { 471 PetscFunctionBegin; 472 if (!*elements) PetscFunctionReturn(0); 473 CHKERRQ(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 CHKERRQ(PetscNew(mesh)); 499 CHKERRQ(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 CHKERRQ(GmshEntitiesDestroy(&(*mesh)->entities)); 510 CHKERRQ(GmshNodesDestroy(&(*mesh)->nodelist)); 511 CHKERRQ(GmshElementsDestroy(&(*mesh)->elements)); 512 CHKERRQ(PetscFree((*mesh)->periodMap)); 513 CHKERRQ(PetscFree((*mesh)->vertexMap)); 514 CHKERRQ(PetscSegBufferDestroy(&(*mesh)->segbuf)); 515 for (r = 0; r < (*mesh)->numRegions; ++r) CHKERRQ(PetscFree((*mesh)->regionNames[r])); 516 CHKERRQ(PetscFree2((*mesh)->regionTags, (*mesh)->regionNames)); 517 CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 531 snum = sscanf(line, "%d", &num); 532 PetscCheckFalse(snum != 1,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 533 CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM)); 539 CHKERRQ(PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE)); 540 if (byteSwap) CHKERRQ(PetscByteSwap(&nid, PETSC_ENUM, 1)); 541 if (byteSwap) CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 565 snum = sscanf(line, "%d", &num); 566 PetscCheckFalse(snum != 1,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 567 CHKERRQ(GmshElementsCreate(num, &elements)); 568 mesh->numElems = num; 569 mesh->elements = elements; 570 for (c = 0; c < num;) { 571 CHKERRQ(PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM)); 572 if (byteSwap) CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, ibuf+off, nint, NULL, PETSC_ENUM)); 587 if (byteSwap) CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG)); 642 if (byteSwap) CHKERRQ(PetscByteSwap(lbuf, PETSC_LONG, 4)); 643 for (i = 0; i < 4; ++i) count[i] = lbuf[i]; 644 CHKERRQ(GmshEntitiesCreate(count, &mesh->entities)); 645 for (dim = 0; dim < 4; ++dim) { 646 for (index = 0; index < count[dim]; ++index) { 647 CHKERRQ(PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM)); 648 if (byteSwap) CHKERRQ(PetscByteSwap(&eid, PETSC_ENUM, 1)); 649 CHKERRQ(GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity)); 650 CHKERRQ(PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE)); 651 if (byteSwap) CHKERRQ(PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6)); 652 CHKERRQ(PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG)); 653 if (byteSwap) CHKERRQ(PetscByteSwap(&num, PETSC_LONG, 1)); 654 CHKERRQ(GmshBufferGet(gmsh, num, sizeof(int), &ibuf)); 655 CHKERRQ(PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM)); 656 if (byteSwap) CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG)); 661 if (byteSwap) CHKERRQ(PetscByteSwap(&num, PETSC_LONG, 1)); 662 CHKERRQ(GmshBufferGet(gmsh, num, sizeof(int), &ibuf)); 663 CHKERRQ(PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM)); 664 if (byteSwap) CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG)); 689 if (byteSwap) CHKERRQ(PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1)); 690 CHKERRQ(PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG)); 691 if (byteSwap) CHKERRQ(PetscByteSwap(&numTotalNodes, PETSC_LONG, 1)); 692 CHKERRQ(GmshNodesCreate(numTotalNodes, &nodes)); 693 mesh->numNodes = numTotalNodes; 694 mesh->nodelist = nodes; 695 for (n = 0, block = 0; block < numEntityBlocks; ++block) { 696 CHKERRQ(PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM)); 697 CHKERRQ(PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG)); 698 if (byteSwap) CHKERRQ(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 CHKERRQ(GmshBufferGet(gmsh, numNodes, nbytes, &cbuf)); 703 CHKERRQ(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()) CHKERRQ(PetscByteSwap(cnid, PETSC_ENUM, 1)); 708 if (!PetscBinaryBigEndian()) CHKERRQ(PetscByteSwap(cxyz, PETSC_DOUBLE, 3)); 709 CHKERRQ(PetscMemcpy(&nid, cnid, sizeof(int))); 710 CHKERRQ(PetscMemcpy(xyz, cxyz, 3*sizeof(double))); 711 if (byteSwap) CHKERRQ(PetscByteSwap(&nid, PETSC_ENUM, 1)); 712 if (byteSwap) CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM)); 720 CHKERRQ(PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE)); 721 if (byteSwap) CHKERRQ(PetscByteSwap(&nid, PETSC_ENUM, 1)); 722 if (byteSwap) CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG)); 753 if (byteSwap) CHKERRQ(PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1)); 754 CHKERRQ(PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG)); 755 if (byteSwap) CHKERRQ(PetscByteSwap(&numTotalElements, PETSC_LONG, 1)); 756 CHKERRQ(GmshElementsCreate(numTotalElements, &elements)); 757 mesh->numElems = numTotalElements; 758 mesh->elements = elements; 759 for (c = 0, block = 0; block < numEntityBlocks; ++block) { 760 CHKERRQ(PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM)); 761 if (byteSwap) CHKERRQ(PetscByteSwap(info, PETSC_ENUM, 3)); 762 eid = info[0]; dim = info[1]; cellType = info[2]; 763 CHKERRQ(GmshEntitiesGet(mesh->entities, dim, eid, &entity)); 764 CHKERRQ(GmshCellTypeCheck(cellType)); 765 numVerts = GmshCellMap[cellType].numVerts; 766 numNodes = GmshCellMap[cellType].numNodes; 767 numTags = entity->numTags; 768 CHKERRQ(PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG)); 769 if (byteSwap) CHKERRQ(PetscByteSwap(&numElements, PETSC_LONG, 1)); 770 CHKERRQ(GmshBufferGet(gmsh, (1+numNodes)*numElements, sizeof(int), &ibuf)); 771 CHKERRQ(PetscViewerRead(viewer, ibuf, (1+numNodes)*numElements, NULL, PETSC_ENUM)); 772 if (byteSwap) CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 803 snum = sscanf(line, "%d", &numPeriodic); 804 PetscCheckFalse(snum != 1,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 805 } else { 806 CHKERRQ(PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM)); 807 if (byteSwap) CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING)); 816 snum = sscanf(line, "%d %d %d", &correspondingDim, &correspondingTag, &primaryTag); 817 PetscCheckFalse(snum != 3,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 818 } else { 819 CHKERRQ(PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM)); 820 if (byteSwap) CHKERRQ(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 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 827 snum = sscanf(line, "%ld", &nNodes); 828 if (snum != 1) { /* discard transformation and try again */ 829 CHKERRQ(PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING)); 830 CHKERRQ(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING)); 831 snum = sscanf(line, "%ld", &nNodes); 832 PetscCheckFalse(snum != 1,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 833 } 834 } else { 835 CHKERRQ(PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG)); 836 if (byteSwap) CHKERRQ(PetscByteSwap(&nNodes, PETSC_LONG, 1)); 837 if (nNodes == -1) { /* discard transformation and try again */ 838 CHKERRQ(PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE)); 839 CHKERRQ(PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG)); 840 if (byteSwap) CHKERRQ(PetscByteSwap(&nNodes, PETSC_LONG, 1)); 841 } 842 } 843 844 for (j = 0; j < nNodes; j++) { 845 if (fileFormat == 22 || !binary) { 846 CHKERRQ(PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING)); 847 snum = sscanf(line, "%d %d", &correspondingNode, &primaryNode); 848 PetscCheckFalse(snum != 2,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 849 } else { 850 CHKERRQ(PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM)); 851 if (byteSwap) CHKERRQ(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 CHKERRQ(GmshReadSize(gmsh, count, 4)); 894 CHKERRQ(GmshEntitiesCreate(count, &mesh->entities)); 895 for (dim = 0; dim < 4; ++dim) { 896 for (index = 0; index < count[dim]; ++index) { 897 CHKERRQ(GmshReadInt(gmsh, &eid, 1)); 898 CHKERRQ(GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity)); 899 CHKERRQ(GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6)); 900 CHKERRQ(GmshReadSize(gmsh, &numTags, 1)); 901 CHKERRQ(GmshBufferGet(gmsh, numTags, sizeof(int), &tags)); 902 CHKERRQ(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 CHKERRQ(GmshReadSize(gmsh, &numTags, 1)); 907 CHKERRQ(GmshBufferGet(gmsh, numTags, sizeof(int), &tags)); 908 CHKERRQ(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 CHKERRQ(GmshReadSize(gmsh, sizes, 4)); 938 numEntityBlocks = sizes[0]; numNodes = sizes[1]; 939 CHKERRQ(GmshNodesCreate(numNodes, &nodes)); 940 mesh->numNodes = numNodes; 941 mesh->nodelist = nodes; 942 for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) { 943 CHKERRQ(GmshReadInt(gmsh, info, 3)); 944 dim = info[0]; eid = info[1]; parametric = info[2]; 945 CHKERRQ(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 CHKERRQ(GmshReadSize(gmsh, &numNodesBlock, 1)); 950 CHKERRQ(GmshReadSize(gmsh, nodes->id+node, numNodesBlock)); 951 CHKERRQ(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 CHKERRQ(GmshReadSize(gmsh, sizes, 4)); 979 numEntityBlocks = sizes[0]; numElements = sizes[1]; 980 CHKERRQ(GmshElementsCreate(numElements, &elements)); 981 mesh->numElems = numElements; 982 mesh->elements = elements; 983 for (c = 0, block = 0; block < numEntityBlocks; ++block) { 984 CHKERRQ(GmshReadInt(gmsh, info, 3)); 985 dim = info[0]; eid = info[1]; cellType = info[2]; 986 CHKERRQ(GmshEntitiesGet(mesh->entities, dim, eid, &entity)); 987 CHKERRQ(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 CHKERRQ(GmshReadSize(gmsh, &numBlockElements, 1)); 993 CHKERRQ(GmshBufferGet(gmsh, (1+numNodes)*numBlockElements, sizeof(PetscInt), &ibuf)); 994 CHKERRQ(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 CHKERRQ(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 CHKERRQ(GmshReadSize(gmsh, &numPeriodicLinks, 1)); 1032 for (link = 0; link < numPeriodicLinks; ++link) { 1033 CHKERRQ(GmshReadInt(gmsh, info, 3)); 1034 CHKERRQ(GmshReadSize(gmsh, &numAffine, 1)); 1035 CHKERRQ(GmshReadDouble(gmsh, dbuf, numAffine)); 1036 CHKERRQ(GmshReadSize(gmsh, &numCorrespondingNodes, 1)); 1037 CHKERRQ(GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags)); 1038 CHKERRQ(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 CHKERRQ(GmshReadString(gmsh, line, 3)); 1064 snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize); 1065 PetscCheckFalse(snum != 3,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line); 1066 PetscCheckFalse(version < 2.2,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version); 1067 PetscCheckFalse((int)version == 3,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version); 1068 PetscCheckFalse(version > 4.1,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version); 1069 PetscCheckFalse(gmsh->binary && !fileType,PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is binary but Gmsh file is ASCII"); 1070 PetscCheckFalse(!gmsh->binary && fileType,PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is ASCII but Gmsh file is binary"); 1071 fileFormat = (int)roundf(version*10); 1072 PetscCheckFalse(fileFormat <= 40 && dataSize != sizeof(double),PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Data size %d is not valid for a Gmsh file", dataSize); 1073 PetscCheckFalse(fileFormat >= 41 && dataSize != sizeof(int) && dataSize != sizeof(PetscInt64),PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Data size %d is not valid for a Gmsh file", dataSize); 1074 gmsh->fileFormat = fileFormat; 1075 gmsh->dataSize = dataSize; 1076 gmsh->byteSwap = PETSC_FALSE; 1077 if (gmsh->binary) { 1078 CHKERRQ(GmshReadInt(gmsh, &checkEndian, 1)); 1079 if (checkEndian != 1) { 1080 CHKERRQ(PetscByteSwap(&checkEndian, PETSC_ENUM, 1)); 1081 PetscCheckFalse(checkEndian != 1,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to detect endianness in Gmsh file header: %s", line); 1082 gmsh->byteSwap = PETSC_TRUE; 1083 } 1084 } 1085 PetscFunctionReturn(0); 1086 } 1087 1088 /* 1089 PhysicalNames 1090 numPhysicalNames(ASCII int) 1091 dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max) 1092 ... 1093 $EndPhysicalNames 1094 */ 1095 static PetscErrorCode GmshReadPhysicalNames(GmshFile *gmsh, GmshMesh *mesh) 1096 { 1097 char line[PETSC_MAX_PATH_LEN], name[128+2], *p, *q; 1098 int snum, region, dim, tag; 1099 1100 PetscFunctionBegin; 1101 CHKERRQ(GmshReadString(gmsh, line, 1)); 1102 snum = sscanf(line, "%d", ®ion); 1103 mesh->numRegions = region; 1104 PetscCheckFalse(snum != 1,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 1105 CHKERRQ(PetscMalloc2(mesh->numRegions, &mesh->regionTags, mesh->numRegions, &mesh->regionNames)); 1106 for (region = 0; region < mesh->numRegions; ++region) { 1107 CHKERRQ(GmshReadString(gmsh, line, 2)); 1108 snum = sscanf(line, "%d %d", &dim, &tag); 1109 PetscCheckFalse(snum != 2,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 1110 CHKERRQ(GmshReadString(gmsh, line, -(PetscInt)sizeof(line))); 1111 CHKERRQ(PetscStrchr(line, '"', &p)); 1112 PetscCheckFalse(!p,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 1113 CHKERRQ(PetscStrrchr(line, '"', &q)); 1114 PetscCheckFalse(q == p,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file"); 1115 CHKERRQ(PetscStrncpy(name, p+1, (size_t)(q-p-1))); 1116 mesh->regionTags[region] = tag; 1117 CHKERRQ(PetscStrallocpy(name, &mesh->regionNames[region])); 1118 } 1119 PetscFunctionReturn(0); 1120 } 1121 1122 static PetscErrorCode GmshReadEntities(GmshFile *gmsh, GmshMesh *mesh) 1123 { 1124 PetscFunctionBegin; 1125 switch (gmsh->fileFormat) { 1126 case 41: CHKERRQ(GmshReadEntities_v41(gmsh, mesh)); break; 1127 default: CHKERRQ(GmshReadEntities_v40(gmsh, mesh)); break; 1128 } 1129 PetscFunctionReturn(0); 1130 } 1131 1132 static PetscErrorCode GmshReadNodes(GmshFile *gmsh, GmshMesh *mesh) 1133 { 1134 PetscFunctionBegin; 1135 switch (gmsh->fileFormat) { 1136 case 41: CHKERRQ(GmshReadNodes_v41(gmsh, mesh)); break; 1137 case 40: CHKERRQ(GmshReadNodes_v40(gmsh, mesh)); break; 1138 default: CHKERRQ(GmshReadNodes_v22(gmsh, mesh)); break; 1139 } 1140 1141 { /* Gmsh v2.2/v4.0 does not provide min/max node tags */ 1142 if (mesh->numNodes > 0 && gmsh->nodeEnd >= gmsh->nodeStart) { 1143 PetscInt tagMin = PETSC_MAX_INT, tagMax = PETSC_MIN_INT, n; 1144 GmshNodes *nodes = mesh->nodelist; 1145 for (n = 0; n < mesh->numNodes; ++n) { 1146 const PetscInt tag = nodes->id[n]; 1147 tagMin = PetscMin(tag, tagMin); 1148 tagMax = PetscMax(tag, tagMax); 1149 } 1150 gmsh->nodeStart = tagMin; 1151 gmsh->nodeEnd = tagMax+1; 1152 } 1153 } 1154 1155 { /* Support for sparse node tags */ 1156 PetscInt n, t; 1157 GmshNodes *nodes = mesh->nodelist; 1158 CHKERRQ(PetscMalloc1(gmsh->nodeEnd - gmsh->nodeStart, &gmsh->nbuf)); 1159 for (t = 0; t < gmsh->nodeEnd - gmsh->nodeStart; ++t) gmsh->nbuf[t] = PETSC_MIN_INT; 1160 gmsh->nodeMap = gmsh->nbuf - gmsh->nodeStart; 1161 for (n = 0; n < mesh->numNodes; ++n) { 1162 const PetscInt tag = nodes->id[n]; 1163 PetscCheckFalse(gmsh->nodeMap[tag] >= 0,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Repeated node tag %D", tag); 1164 gmsh->nodeMap[tag] = n; 1165 } 1166 } 1167 PetscFunctionReturn(0); 1168 } 1169 1170 static PetscErrorCode GmshReadElements(GmshFile *gmsh, GmshMesh *mesh) 1171 { 1172 PetscFunctionBegin; 1173 switch (gmsh->fileFormat) { 1174 case 41: CHKERRQ(GmshReadElements_v41(gmsh, mesh)); break; 1175 case 40: CHKERRQ(GmshReadElements_v40(gmsh, mesh)); break; 1176 default: CHKERRQ(GmshReadElements_v22(gmsh, mesh)); break; 1177 } 1178 1179 { /* Reorder elements by codimension and polytope type */ 1180 PetscInt ne = mesh->numElems; 1181 GmshElement *elements = mesh->elements; 1182 PetscInt keymap[GMSH_NUM_POLYTOPES], nk = 0; 1183 PetscInt offset[GMSH_NUM_POLYTOPES+1], e, k; 1184 1185 for (k = 0; k < GMSH_NUM_POLYTOPES; ++k) keymap[k] = PETSC_MIN_INT; 1186 CHKERRQ(PetscMemzero(offset,sizeof(offset))); 1187 1188 keymap[GMSH_TET] = nk++; 1189 keymap[GMSH_HEX] = nk++; 1190 keymap[GMSH_PRI] = nk++; 1191 keymap[GMSH_PYR] = nk++; 1192 keymap[GMSH_TRI] = nk++; 1193 keymap[GMSH_QUA] = nk++; 1194 keymap[GMSH_SEG] = nk++; 1195 keymap[GMSH_VTX] = nk++; 1196 1197 CHKERRQ(GmshElementsCreate(mesh->numElems, &mesh->elements)); 1198 #define key(eid) keymap[GmshCellMap[elements[(eid)].cellType].polytope] 1199 for (e = 0; e < ne; ++e) offset[1+key(e)]++; 1200 for (k = 1; k < nk; ++k) offset[k] += offset[k-1]; 1201 for (e = 0; e < ne; ++e) mesh->elements[offset[key(e)]++] = elements[e]; 1202 #undef key 1203 CHKERRQ(GmshElementsDestroy(&elements)); 1204 } 1205 1206 { /* Mesh dimension and order */ 1207 GmshElement *elem = mesh->numElems ? mesh->elements : NULL; 1208 mesh->dim = elem ? GmshCellMap[elem->cellType].dim : 0; 1209 mesh->order = elem ? GmshCellMap[elem->cellType].order : 0; 1210 } 1211 1212 { 1213 PetscBT vtx; 1214 PetscInt dim = mesh->dim, e, n, v; 1215 1216 CHKERRQ(PetscBTCreate(mesh->numNodes, &vtx)); 1217 1218 /* Compute number of cells and set of vertices */ 1219 mesh->numCells = 0; 1220 for (e = 0; e < mesh->numElems; ++e) { 1221 GmshElement *elem = mesh->elements + e; 1222 if (elem->dim == dim && dim > 0) mesh->numCells++; 1223 for (v = 0; v < elem->numVerts; v++) { 1224 CHKERRQ(PetscBTSet(vtx, elem->nodes[v])); 1225 } 1226 } 1227 1228 /* Compute numbering for vertices */ 1229 mesh->numVerts = 0; 1230 CHKERRQ(PetscMalloc1(mesh->numNodes, &mesh->vertexMap)); 1231 for (n = 0; n < mesh->numNodes; ++n) 1232 mesh->vertexMap[n] = PetscBTLookup(vtx, n) ? mesh->numVerts++ : PETSC_MIN_INT; 1233 1234 CHKERRQ(PetscBTDestroy(&vtx)); 1235 } 1236 PetscFunctionReturn(0); 1237 } 1238 1239 static PetscErrorCode GmshReadPeriodic(GmshFile *gmsh, GmshMesh *mesh) 1240 { 1241 PetscInt n; 1242 1243 PetscFunctionBegin; 1244 CHKERRQ(PetscMalloc1(mesh->numNodes, &mesh->periodMap)); 1245 for (n = 0; n < mesh->numNodes; ++n) mesh->periodMap[n] = n; 1246 switch (gmsh->fileFormat) { 1247 case 41: CHKERRQ(GmshReadPeriodic_v41(gmsh, mesh->periodMap)); break; 1248 default: CHKERRQ(GmshReadPeriodic_v40(gmsh, mesh->periodMap)); break; 1249 } 1250 1251 /* Find canonical primary nodes */ 1252 for (n = 0; n < mesh->numNodes; ++n) 1253 while (mesh->periodMap[n] != mesh->periodMap[mesh->periodMap[n]]) 1254 mesh->periodMap[n] = mesh->periodMap[mesh->periodMap[n]]; 1255 1256 /* Renumber vertices (filter out correspondings) */ 1257 mesh->numVerts = 0; 1258 for (n = 0; n < mesh->numNodes; ++n) 1259 if (mesh->vertexMap[n] >= 0) /* is vertex */ 1260 if (mesh->periodMap[n] == n) /* is primary */ 1261 mesh->vertexMap[n] = mesh->numVerts++; 1262 for (n = 0; n < mesh->numNodes; ++n) 1263 if (mesh->vertexMap[n] >= 0) /* is vertex */ 1264 if (mesh->periodMap[n] != n) /* is corresponding */ 1265 mesh->vertexMap[n] = mesh->vertexMap[mesh->periodMap[n]]; 1266 PetscFunctionReturn(0); 1267 } 1268 1269 #define DM_POLYTOPE_VERTEX DM_POLYTOPE_POINT 1270 #define DM_POLYTOPE_PYRAMID DM_POLYTOPE_UNKNOWN 1271 static const DMPolytopeType DMPolytopeMap[] = { 1272 /* GMSH_VTX */ DM_POLYTOPE_VERTEX, 1273 /* GMSH_SEG */ DM_POLYTOPE_SEGMENT, 1274 /* GMSH_TRI */ DM_POLYTOPE_TRIANGLE, 1275 /* GMSH_QUA */ DM_POLYTOPE_QUADRILATERAL, 1276 /* GMSH_TET */ DM_POLYTOPE_TETRAHEDRON, 1277 /* GMSH_HEX */ DM_POLYTOPE_HEXAHEDRON, 1278 /* GMSH_PRI */ DM_POLYTOPE_TRI_PRISM, 1279 /* GMSH_PYR */ DM_POLYTOPE_PYRAMID, 1280 DM_POLYTOPE_UNKNOWN 1281 }; 1282 1283 static inline DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt cellType) 1284 { 1285 return DMPolytopeMap[GmshCellMap[cellType].polytope]; 1286 } 1287 1288 static PetscErrorCode GmshCreateFE(MPI_Comm comm, const char prefix[], PetscBool isSimplex, PetscBool continuity, PetscDTNodeType nodeType, PetscInt dim, PetscInt Nc, PetscInt k, PetscFE *fem) 1289 { 1290 DM K; 1291 PetscSpace P; 1292 PetscDualSpace Q; 1293 PetscQuadrature q, fq; 1294 PetscBool isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 1295 PetscBool endpoint = PETSC_TRUE; 1296 char name[32]; 1297 1298 PetscFunctionBegin; 1299 /* Create space */ 1300 CHKERRQ(PetscSpaceCreate(comm, &P)); 1301 CHKERRQ(PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL)); 1302 CHKERRQ(PetscSpacePolynomialSetTensor(P, isTensor)); 1303 CHKERRQ(PetscSpaceSetNumComponents(P, Nc)); 1304 CHKERRQ(PetscSpaceSetNumVariables(P, dim)); 1305 CHKERRQ(PetscSpaceSetDegree(P, k, PETSC_DETERMINE)); 1306 if (prefix) { 1307 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject) P, prefix)); 1308 CHKERRQ(PetscSpaceSetFromOptions(P)); 1309 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject) P, NULL)); 1310 CHKERRQ(PetscSpaceGetDegree(P, &k, NULL)); 1311 } 1312 CHKERRQ(PetscSpaceSetUp(P)); 1313 /* Create dual space */ 1314 CHKERRQ(PetscDualSpaceCreate(comm, &Q)); 1315 CHKERRQ(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE)); 1316 CHKERRQ(PetscDualSpaceLagrangeSetTensor(Q, isTensor)); 1317 CHKERRQ(PetscDualSpaceLagrangeSetContinuity(Q, continuity)); 1318 CHKERRQ(PetscDualSpaceLagrangeSetNodeType(Q, nodeType, endpoint, 0)); 1319 CHKERRQ(PetscDualSpaceSetNumComponents(Q, Nc)); 1320 CHKERRQ(PetscDualSpaceSetOrder(Q, k)); 1321 CHKERRQ(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K)); 1322 CHKERRQ(PetscDualSpaceSetDM(Q, K)); 1323 CHKERRQ(DMDestroy(&K)); 1324 if (prefix) { 1325 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject) Q, prefix)); 1326 CHKERRQ(PetscDualSpaceSetFromOptions(Q)); 1327 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject) Q, NULL)); 1328 } 1329 CHKERRQ(PetscDualSpaceSetUp(Q)); 1330 /* Create quadrature */ 1331 if (isSimplex) { 1332 CHKERRQ(PetscDTStroudConicalQuadrature(dim, 1, k+1, -1, +1, &q)); 1333 CHKERRQ(PetscDTStroudConicalQuadrature(dim-1, 1, k+1, -1, +1, &fq)); 1334 } else { 1335 CHKERRQ(PetscDTGaussTensorQuadrature(dim, 1, k+1, -1, +1, &q)); 1336 CHKERRQ(PetscDTGaussTensorQuadrature(dim-1, 1, k+1, -1, +1, &fq)); 1337 } 1338 /* Create finite element */ 1339 CHKERRQ(PetscFECreate(comm, fem)); 1340 CHKERRQ(PetscFESetType(*fem, PETSCFEBASIC)); 1341 CHKERRQ(PetscFESetNumComponents(*fem, Nc)); 1342 CHKERRQ(PetscFESetBasisSpace(*fem, P)); 1343 CHKERRQ(PetscFESetDualSpace(*fem, Q)); 1344 CHKERRQ(PetscFESetQuadrature(*fem, q)); 1345 CHKERRQ(PetscFESetFaceQuadrature(*fem, fq)); 1346 if (prefix) { 1347 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix)); 1348 CHKERRQ(PetscFESetFromOptions(*fem)); 1349 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject) *fem, NULL)); 1350 } 1351 CHKERRQ(PetscFESetUp(*fem)); 1352 /* Cleanup */ 1353 CHKERRQ(PetscSpaceDestroy(&P)); 1354 CHKERRQ(PetscDualSpaceDestroy(&Q)); 1355 CHKERRQ(PetscQuadratureDestroy(&q)); 1356 CHKERRQ(PetscQuadratureDestroy(&fq)); 1357 /* Set finite element name */ 1358 CHKERRQ(PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k)); 1359 CHKERRQ(PetscFESetName(*fem, name)); 1360 PetscFunctionReturn(0); 1361 } 1362 1363 /*@C 1364 DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file 1365 1366 + comm - The MPI communicator 1367 . filename - Name of the Gmsh file 1368 - interpolate - Create faces and edges in the mesh 1369 1370 Output Parameter: 1371 . dm - The DM object representing the mesh 1372 1373 Level: beginner 1374 1375 .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate() 1376 @*/ 1377 PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 1378 { 1379 PetscViewer viewer; 1380 PetscMPIInt rank; 1381 int fileType; 1382 PetscViewerType vtype; 1383 1384 PetscFunctionBegin; 1385 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 1386 1387 /* Determine Gmsh file type (ASCII or binary) from file header */ 1388 if (rank == 0) { 1389 GmshFile gmsh[1]; 1390 char line[PETSC_MAX_PATH_LEN]; 1391 int snum; 1392 float version; 1393 1394 CHKERRQ(PetscArrayzero(gmsh,1)); 1395 CHKERRQ(PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer)); 1396 CHKERRQ(PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII)); 1397 CHKERRQ(PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ)); 1398 CHKERRQ(PetscViewerFileSetName(gmsh->viewer, filename)); 1399 /* Read only the first two lines of the Gmsh file */ 1400 CHKERRQ(GmshReadSection(gmsh, line)); 1401 CHKERRQ(GmshExpect(gmsh, "$MeshFormat", line)); 1402 CHKERRQ(GmshReadString(gmsh, line, 2)); 1403 snum = sscanf(line, "%f %d", &version, &fileType); 1404 PetscCheckFalse(snum != 2,PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line); 1405 PetscCheckFalse(version < 2.2,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version); 1406 PetscCheckFalse((int)version == 3,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version); 1407 PetscCheckFalse(version > 4.1,PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version); 1408 CHKERRQ(PetscViewerDestroy(&gmsh->viewer)); 1409 } 1410 CHKERRMPI(MPI_Bcast(&fileType, 1, MPI_INT, 0, comm)); 1411 vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY; 1412 1413 /* Create appropriate viewer and build plex */ 1414 CHKERRQ(PetscViewerCreate(comm, &viewer)); 1415 CHKERRQ(PetscViewerSetType(viewer, vtype)); 1416 CHKERRQ(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 1417 CHKERRQ(PetscViewerFileSetName(viewer, filename)); 1418 CHKERRQ(DMPlexCreateGmsh(comm, viewer, interpolate, dm)); 1419 CHKERRQ(PetscViewerDestroy(&viewer)); 1420 PetscFunctionReturn(0); 1421 } 1422 1423 /*@ 1424 DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer 1425 1426 Collective 1427 1428 Input Parameters: 1429 + comm - The MPI communicator 1430 . viewer - The Viewer associated with a Gmsh file 1431 - interpolate - Create faces and edges in the mesh 1432 1433 Output Parameter: 1434 . dm - The DM object representing the mesh 1435 1436 Note: http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format 1437 1438 Level: beginner 1439 1440 .seealso: DMPLEX, DMCreate() 1441 @*/ 1442 PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm) 1443 { 1444 GmshMesh *mesh = NULL; 1445 PetscViewer parentviewer = NULL; 1446 PetscBT periodicVerts = NULL; 1447 PetscBT periodicCells = NULL; 1448 DM cdm; 1449 PetscSection coordSection; 1450 Vec coordinates; 1451 DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL, marker = NULL, *regionSets; 1452 PetscInt dim = 0, coordDim = -1, order = 0; 1453 PetscInt numNodes = 0, numElems = 0, numVerts = 0, numCells = 0; 1454 PetscInt cell, cone[8], e, n, v, d; 1455 PetscBool binary, usemarker = PETSC_FALSE, useregions = PETSC_FALSE, markvertices = PETSC_FALSE; 1456 PetscBool hybrid = interpolate, periodic = PETSC_TRUE; 1457 PetscBool highOrder = PETSC_TRUE, highOrderSet, project = PETSC_FALSE; 1458 PetscBool isSimplex = PETSC_FALSE, isHybrid = PETSC_FALSE, hasTetra = PETSC_FALSE; 1459 PetscMPIInt rank; 1460 PetscErrorCode ierr; 1461 1462 PetscFunctionBegin; 1463 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 1464 ierr = PetscObjectOptionsBegin((PetscObject)viewer);CHKERRQ(ierr); 1465 CHKERRQ(PetscOptionsHead(PetscOptionsObject,"DMPlex Gmsh options")); 1466 CHKERRQ(PetscOptionsBool("-dm_plex_gmsh_hybrid", "Generate hybrid cell bounds", "DMPlexCreateGmsh", hybrid, &hybrid, NULL)); 1467 CHKERRQ(PetscOptionsBool("-dm_plex_gmsh_periodic","Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL)); 1468 CHKERRQ(PetscOptionsBool("-dm_plex_gmsh_highorder","Generate high-order coordinates", "DMPlexCreateGmsh", highOrder, &highOrder, &highOrderSet)); 1469 CHKERRQ(PetscOptionsBool("-dm_plex_gmsh_project", "Project high-order coordinates to a different space", "DMPlexCreateGmsh", project, &project, NULL)); 1470 CHKERRQ(PetscOptionsBool("-dm_plex_gmsh_use_marker", "Generate marker label", "DMPlexCreateGmsh", usemarker, &usemarker, NULL)); 1471 CHKERRQ(PetscOptionsBool("-dm_plex_gmsh_use_regions", "Generate labels with region names", "DMPlexCreateGmsh", useregions, &useregions, NULL)); 1472 CHKERRQ(PetscOptionsBool("-dm_plex_gmsh_mark_vertices", "Add vertices to generated labels", "DMPlexCreateGmsh", markvertices, &markvertices, NULL)); 1473 CHKERRQ(PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE)); 1474 CHKERRQ(PetscOptionsTail()); 1475 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1476 1477 CHKERRQ(GmshCellInfoSetUp()); 1478 1479 CHKERRQ(DMCreate(comm, dm)); 1480 CHKERRQ(DMSetType(*dm, DMPLEX)); 1481 CHKERRQ(PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,NULL,NULL,NULL)); 1482 1483 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary)); 1484 1485 /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */ 1486 if (binary) { 1487 parentviewer = viewer; 1488 CHKERRQ(PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer)); 1489 } 1490 1491 if (rank == 0) { 1492 GmshFile gmsh[1]; 1493 char line[PETSC_MAX_PATH_LEN]; 1494 PetscBool match; 1495 1496 CHKERRQ(PetscArrayzero(gmsh,1)); 1497 gmsh->viewer = viewer; 1498 gmsh->binary = binary; 1499 1500 CHKERRQ(GmshMeshCreate(&mesh)); 1501 1502 /* Read mesh format */ 1503 CHKERRQ(GmshReadSection(gmsh, line)); 1504 CHKERRQ(GmshExpect(gmsh, "$MeshFormat", line)); 1505 CHKERRQ(GmshReadMeshFormat(gmsh)); 1506 CHKERRQ(GmshReadEndSection(gmsh, "$EndMeshFormat", line)); 1507 1508 /* OPTIONAL Read physical names */ 1509 CHKERRQ(GmshReadSection(gmsh, line)); 1510 CHKERRQ(GmshMatch(gmsh, "$PhysicalNames", line, &match)); 1511 if (match) { 1512 CHKERRQ(GmshExpect(gmsh, "$PhysicalNames", line)); 1513 CHKERRQ(GmshReadPhysicalNames(gmsh, mesh)); 1514 CHKERRQ(GmshReadEndSection(gmsh, "$EndPhysicalNames", line)); 1515 /* Initial read for entity section */ 1516 CHKERRQ(GmshReadSection(gmsh, line)); 1517 } 1518 1519 /* Read entities */ 1520 if (gmsh->fileFormat >= 40) { 1521 CHKERRQ(GmshExpect(gmsh, "$Entities", line)); 1522 CHKERRQ(GmshReadEntities(gmsh, mesh)); 1523 CHKERRQ(GmshReadEndSection(gmsh, "$EndEntities", line)); 1524 /* Initial read for nodes section */ 1525 CHKERRQ(GmshReadSection(gmsh, line)); 1526 } 1527 1528 /* Read nodes */ 1529 CHKERRQ(GmshExpect(gmsh, "$Nodes", line)); 1530 CHKERRQ(GmshReadNodes(gmsh, mesh)); 1531 CHKERRQ(GmshReadEndSection(gmsh, "$EndNodes", line)); 1532 1533 /* Read elements */ 1534 CHKERRQ(GmshReadSection(gmsh, line)); 1535 CHKERRQ(GmshExpect(gmsh, "$Elements", line)); 1536 CHKERRQ(GmshReadElements(gmsh, mesh)); 1537 CHKERRQ(GmshReadEndSection(gmsh, "$EndElements", line)); 1538 1539 /* Read periodic section (OPTIONAL) */ 1540 if (periodic) { 1541 CHKERRQ(GmshReadSection(gmsh, line)); 1542 CHKERRQ(GmshMatch(gmsh, "$Periodic", line, &periodic)); 1543 } 1544 if (periodic) { 1545 CHKERRQ(GmshExpect(gmsh, "$Periodic", line)); 1546 CHKERRQ(GmshReadPeriodic(gmsh, mesh)); 1547 CHKERRQ(GmshReadEndSection(gmsh, "$EndPeriodic", line)); 1548 } 1549 1550 CHKERRQ(PetscFree(gmsh->wbuf)); 1551 CHKERRQ(PetscFree(gmsh->sbuf)); 1552 CHKERRQ(PetscFree(gmsh->nbuf)); 1553 1554 dim = mesh->dim; 1555 order = mesh->order; 1556 numNodes = mesh->numNodes; 1557 numElems = mesh->numElems; 1558 numVerts = mesh->numVerts; 1559 numCells = mesh->numCells; 1560 1561 { 1562 GmshElement *elemA = mesh->numCells > 0 ? mesh->elements : NULL; 1563 GmshElement *elemB = elemA ? elemA + mesh->numCells - 1 : NULL; 1564 int ptA = elemA ? GmshCellMap[elemA->cellType].polytope : -1; 1565 int ptB = elemB ? GmshCellMap[elemB->cellType].polytope : -1; 1566 isSimplex = (ptA == GMSH_QUA || ptA == GMSH_HEX) ? PETSC_FALSE : PETSC_TRUE; 1567 isHybrid = (ptA == ptB) ? PETSC_FALSE : PETSC_TRUE; 1568 hasTetra = (ptA == GMSH_TET) ? PETSC_TRUE : PETSC_FALSE; 1569 } 1570 } 1571 1572 if (parentviewer) { 1573 CHKERRQ(PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer)); 1574 } 1575 1576 { 1577 int buf[6]; 1578 1579 buf[0] = (int)dim; 1580 buf[1] = (int)order; 1581 buf[2] = periodic; 1582 buf[3] = isSimplex; 1583 buf[4] = isHybrid; 1584 buf[5] = hasTetra; 1585 1586 CHKERRMPI(MPI_Bcast(buf, 6, MPI_INT, 0, comm)); 1587 1588 dim = buf[0]; 1589 order = buf[1]; 1590 periodic = buf[2] ? PETSC_TRUE : PETSC_FALSE; 1591 isSimplex = buf[3] ? PETSC_TRUE : PETSC_FALSE; 1592 isHybrid = buf[4] ? PETSC_TRUE : PETSC_FALSE; 1593 hasTetra = buf[5] ? PETSC_TRUE : PETSC_FALSE; 1594 } 1595 1596 if (!highOrderSet) highOrder = (order > 1) ? PETSC_TRUE : PETSC_FALSE; 1597 PetscCheckFalse(highOrder && isHybrid,comm, PETSC_ERR_SUP, "No support for discretization on hybrid meshes yet"); 1598 1599 /* We do not want this label automatically computed, instead we fill it here */ 1600 CHKERRQ(DMCreateLabel(*dm, "celltype")); 1601 1602 /* Allocate the cell-vertex mesh */ 1603 CHKERRQ(DMPlexSetChart(*dm, 0, numCells+numVerts)); 1604 for (cell = 0; cell < numCells; ++cell) { 1605 GmshElement *elem = mesh->elements + cell; 1606 DMPolytopeType ctype = DMPolytopeTypeFromGmsh(elem->cellType); 1607 if (hybrid && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR; 1608 CHKERRQ(DMPlexSetConeSize(*dm, cell, elem->numVerts)); 1609 CHKERRQ(DMPlexSetCellType(*dm, cell, ctype)); 1610 } 1611 for (v = numCells; v < numCells+numVerts; ++v) { 1612 CHKERRQ(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT)); 1613 } 1614 CHKERRQ(DMSetUp(*dm)); 1615 1616 /* Add cell-vertex connections */ 1617 for (cell = 0; cell < numCells; ++cell) { 1618 GmshElement *elem = mesh->elements + cell; 1619 for (v = 0; v < elem->numVerts; ++v) { 1620 const PetscInt nn = elem->nodes[v]; 1621 const PetscInt vv = mesh->vertexMap[nn]; 1622 cone[v] = numCells + vv; 1623 } 1624 CHKERRQ(DMPlexReorderCell(*dm, cell, cone)); 1625 CHKERRQ(DMPlexSetCone(*dm, cell, cone)); 1626 } 1627 1628 CHKERRQ(DMSetDimension(*dm, dim)); 1629 CHKERRQ(DMPlexSymmetrize(*dm)); 1630 CHKERRQ(DMPlexStratify(*dm)); 1631 if (interpolate) { 1632 DM idm; 1633 1634 CHKERRQ(DMPlexInterpolate(*dm, &idm)); 1635 CHKERRQ(DMDestroy(dm)); 1636 *dm = idm; 1637 } 1638 1639 /* Create the label "marker" over the whole boundary */ 1640 PetscCheckFalse(usemarker && !interpolate && dim > 1,comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation"); 1641 if (rank == 0 && usemarker) { 1642 PetscInt f, fStart, fEnd; 1643 1644 CHKERRQ(DMCreateLabel(*dm, "marker")); 1645 CHKERRQ(DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd)); 1646 for (f = fStart; f < fEnd; ++f) { 1647 PetscInt suppSize; 1648 1649 CHKERRQ(DMPlexGetSupportSize(*dm, f, &suppSize)); 1650 if (suppSize == 1) { 1651 PetscInt *cone = NULL, coneSize, p; 1652 1653 CHKERRQ(DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone)); 1654 for (p = 0; p < coneSize; p += 2) { 1655 CHKERRQ(DMSetLabelValue_Fast(*dm, &marker, "marker", cone[p], 1)); 1656 } 1657 CHKERRQ(DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone)); 1658 } 1659 } 1660 } 1661 1662 if (rank == 0) { 1663 const PetscInt Nr = useregions ? mesh->numRegions : 0; 1664 PetscInt vStart, vEnd; 1665 1666 CHKERRQ(PetscCalloc1(Nr, ®ionSets)); 1667 CHKERRQ(DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd)); 1668 for (cell = 0, e = 0; e < numElems; ++e) { 1669 GmshElement *elem = mesh->elements + e; 1670 1671 /* Create cell sets */ 1672 if (elem->dim == dim && dim > 0) { 1673 if (elem->numTags > 0) { 1674 const PetscInt tag = elem->tags[0]; 1675 PetscInt r; 1676 1677 if (!Nr) CHKERRQ(DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", cell, tag)); 1678 for (r = 0; r < Nr; ++r) { 1679 if (mesh->regionTags[r] == tag) CHKERRQ(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], cell, tag)); 1680 } 1681 } 1682 cell++; 1683 } 1684 1685 /* Create face sets */ 1686 if (interpolate && elem->dim == dim-1) { 1687 PetscInt joinSize; 1688 const PetscInt *join = NULL; 1689 const PetscInt tag = elem->tags[0]; 1690 PetscInt r; 1691 1692 /* Find the relevant facet with vertex joins */ 1693 for (v = 0; v < elem->numVerts; ++v) { 1694 const PetscInt nn = elem->nodes[v]; 1695 const PetscInt vv = mesh->vertexMap[nn]; 1696 cone[v] = vStart + vv; 1697 } 1698 CHKERRQ(DMPlexGetFullJoin(*dm, elem->numVerts, cone, &joinSize, &join)); 1699 PetscCheckFalse(joinSize != 1,PETSC_COMM_SELF, PETSC_ERR_SUP, "Could not determine Plex facet for Gmsh element %D (Plex cell %D)", elem->id, e); 1700 if (!Nr) CHKERRQ(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], tag)); 1701 for (r = 0; r < Nr; ++r) { 1702 if (mesh->regionTags[r] == tag) CHKERRQ(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], join[0], tag)); 1703 } 1704 CHKERRQ(DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join)); 1705 } 1706 1707 /* Create vertex sets */ 1708 if (elem->dim == 0) { 1709 if (elem->numTags > 0) { 1710 const PetscInt nn = elem->nodes[0]; 1711 const PetscInt vv = mesh->vertexMap[nn]; 1712 const PetscInt tag = elem->tags[0]; 1713 PetscInt r; 1714 1715 if (!Nr) CHKERRQ(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag)); 1716 for (r = 0; r < Nr; ++r) { 1717 if (mesh->regionTags[r] == tag) CHKERRQ(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag)); 1718 } 1719 } 1720 } 1721 } 1722 if (markvertices) { 1723 for (v = 0; v < numNodes; ++v) { 1724 const PetscInt vv = mesh->vertexMap[v]; 1725 const PetscInt tag = mesh->nodelist->tag[v]; 1726 PetscInt r; 1727 1728 if (tag != -1) { 1729 if (!Nr) CHKERRQ(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag)); 1730 for (r = 0; r < Nr; ++r) { 1731 if (mesh->regionTags[r] == tag) CHKERRQ(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag)); 1732 } 1733 } 1734 } 1735 } 1736 CHKERRQ(PetscFree(regionSets)); 1737 } 1738 1739 { /* Create Cell/Face/Vertex Sets labels at all processes */ 1740 enum {n = 4}; 1741 PetscBool flag[n]; 1742 1743 flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE; 1744 flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE; 1745 flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE; 1746 flag[3] = marker ? PETSC_TRUE : PETSC_FALSE; 1747 CHKERRMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm)); 1748 if (flag[0]) CHKERRQ(DMCreateLabel(*dm, "Cell Sets")); 1749 if (flag[1]) CHKERRQ(DMCreateLabel(*dm, "Face Sets")); 1750 if (flag[2]) CHKERRQ(DMCreateLabel(*dm, "Vertex Sets")); 1751 if (flag[3]) CHKERRQ(DMCreateLabel(*dm, "marker")); 1752 } 1753 1754 if (periodic) { 1755 CHKERRQ(PetscBTCreate(numVerts, &periodicVerts)); 1756 for (n = 0; n < numNodes; ++n) { 1757 if (mesh->vertexMap[n] >= 0) { 1758 if (PetscUnlikely(mesh->periodMap[n] != n)) { 1759 PetscInt m = mesh->periodMap[n]; 1760 CHKERRQ(PetscBTSet(periodicVerts, mesh->vertexMap[n])); 1761 CHKERRQ(PetscBTSet(periodicVerts, mesh->vertexMap[m])); 1762 } 1763 } 1764 } 1765 CHKERRQ(PetscBTCreate(numCells, &periodicCells)); 1766 for (cell = 0; cell < numCells; ++cell) { 1767 GmshElement *elem = mesh->elements + cell; 1768 for (v = 0; v < elem->numVerts; ++v) { 1769 PetscInt nn = elem->nodes[v]; 1770 PetscInt vv = mesh->vertexMap[nn]; 1771 if (PetscUnlikely(PetscBTLookup(periodicVerts, vv))) { 1772 CHKERRQ(PetscBTSet(periodicCells, cell)); break; 1773 } 1774 } 1775 } 1776 } 1777 1778 /* Setup coordinate DM */ 1779 if (coordDim < 0) coordDim = dim; 1780 CHKERRQ(DMSetCoordinateDim(*dm, coordDim)); 1781 CHKERRQ(DMGetCoordinateDM(*dm, &cdm)); 1782 if (highOrder) { 1783 PetscFE fe; 1784 PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE; 1785 PetscDTNodeType nodeType = PETSCDTNODES_EQUISPACED; 1786 1787 if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */ 1788 1789 CHKERRQ(GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe)); 1790 CHKERRQ(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view")); 1791 CHKERRQ(DMSetField(cdm, 0, NULL, (PetscObject) fe)); 1792 CHKERRQ(PetscFEDestroy(&fe)); 1793 CHKERRQ(DMCreateDS(cdm)); 1794 } 1795 1796 /* Create coordinates */ 1797 if (highOrder) { 1798 1799 PetscInt maxDof = GmshNumNodes_HEX(order)*coordDim; 1800 double *coords = mesh ? mesh->nodelist->xyz : NULL; 1801 PetscSection section; 1802 PetscScalar *cellCoords; 1803 1804 CHKERRQ(DMSetLocalSection(cdm, NULL)); 1805 CHKERRQ(DMGetLocalSection(cdm, &coordSection)); 1806 CHKERRQ(PetscSectionClone(coordSection, §ion)); 1807 CHKERRQ(DMPlexSetClosurePermutationTensor(cdm, 0, section)); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */ 1808 1809 CHKERRQ(DMCreateLocalVector(cdm, &coordinates)); 1810 CHKERRQ(PetscMalloc1(maxDof, &cellCoords)); 1811 for (cell = 0; cell < numCells; ++cell) { 1812 GmshElement *elem = mesh->elements + cell; 1813 const int *lexorder = GmshCellMap[elem->cellType].lexorder(); 1814 int s = 0; 1815 for (n = 0; n < elem->numNodes; ++n) { 1816 while (lexorder[n+s] < 0) ++s; 1817 const PetscInt node = elem->nodes[lexorder[n+s]]; 1818 for (d = 0; d < coordDim; ++d) cellCoords[(n+s)*coordDim+d] = (PetscReal) coords[node*3+d]; 1819 } 1820 if (s) { 1821 /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */ 1822 PetscReal quaCenterWeights[9] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25}; 1823 /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */ 1824 PetscReal hexBottomWeights[27] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25, 1825 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1826 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; 1827 PetscReal hexFrontWeights[27] = {-0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1828 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1829 -0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; 1830 PetscReal hexLeftWeights[27] = {-0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 1831 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 1832 -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0}; 1833 PetscReal hexRightWeights[27] = { 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 1834 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1835 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25}; 1836 PetscReal hexBackWeights[27] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25, 1837 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 1838 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25}; 1839 PetscReal hexTopWeights[27] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1840 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1841 -0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25}; 1842 PetscReal hexCenterWeights[27] = {-0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25, 1843 0.25, 0.0, 0.25, 0.0, 0.0, 0.0, 0.25, 0.0, 0.25, 1844 -0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25}; 1845 PetscReal *sdWeights2[9] = {NULL, NULL, NULL, NULL, quaCenterWeights, NULL, NULL, NULL, NULL}; 1846 PetscReal *sdWeights3[27] = {NULL, NULL, NULL, NULL, hexBottomWeights, NULL, NULL, NULL, NULL, 1847 NULL, hexFrontWeights, NULL, hexLeftWeights, hexCenterWeights, hexRightWeights, NULL, hexBackWeights, NULL, 1848 NULL, NULL, NULL, NULL, hexTopWeights, NULL, NULL, NULL, NULL}; 1849 PetscReal **sdWeights[4] = {NULL, NULL, sdWeights2, sdWeights3}; 1850 1851 /* Missing entries in serendipity cell, only works for 8-node quad and 20-node hex */ 1852 for (n = 0; n < elem->numNodes+s; ++n) { 1853 if (lexorder[n] >= 0) continue; 1854 for (d = 0; d < coordDim; ++d) cellCoords[n*coordDim+d] = 0.0; 1855 for (int bn = 0; bn < elem->numNodes+s; ++bn) { 1856 if (lexorder[bn] < 0) continue; 1857 const PetscReal *weights = sdWeights[coordDim][n]; 1858 const PetscInt bnode = elem->nodes[lexorder[bn]]; 1859 for (d = 0; d < coordDim; ++d) cellCoords[n*coordDim+d] += weights[bn] * (PetscReal) coords[bnode*3+d]; 1860 } 1861 } 1862 } 1863 CHKERRQ(DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES)); 1864 } 1865 CHKERRQ(PetscSectionDestroy(§ion)); 1866 CHKERRQ(PetscFree(cellCoords)); 1867 1868 } else { 1869 1870 PetscInt *nodeMap; 1871 double *coords = mesh ? mesh->nodelist->xyz : NULL; 1872 PetscScalar *pointCoords; 1873 1874 CHKERRQ(DMGetLocalSection(cdm, &coordSection)); 1875 CHKERRQ(PetscSectionSetNumFields(coordSection, 1)); 1876 CHKERRQ(PetscSectionSetFieldComponents(coordSection, 0, coordDim)); 1877 if (periodic) { /* we need to localize coordinates on cells */ 1878 CHKERRQ(PetscSectionSetChart(coordSection, 0, numCells+numVerts)); 1879 } else { 1880 CHKERRQ(PetscSectionSetChart(coordSection, numCells, numCells+numVerts)); 1881 } 1882 if (periodic) { 1883 for (cell = 0; cell < numCells; ++cell) { 1884 if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) { 1885 GmshElement *elem = mesh->elements + cell; 1886 PetscInt dof = elem->numVerts * coordDim; 1887 CHKERRQ(PetscSectionSetDof(coordSection, cell, dof)); 1888 CHKERRQ(PetscSectionSetFieldDof(coordSection, cell, 0, dof)); 1889 } 1890 } 1891 } 1892 for (v = numCells; v < numCells+numVerts; ++v) { 1893 CHKERRQ(PetscSectionSetDof(coordSection, v, coordDim)); 1894 CHKERRQ(PetscSectionSetFieldDof(coordSection, v, 0, coordDim)); 1895 } 1896 CHKERRQ(PetscSectionSetUp(coordSection)); 1897 1898 CHKERRQ(DMCreateLocalVector(cdm, &coordinates)); 1899 CHKERRQ(VecGetArray(coordinates, &pointCoords)); 1900 if (periodic) { 1901 for (cell = 0; cell < numCells; ++cell) { 1902 if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) { 1903 GmshElement *elem = mesh->elements + cell; 1904 PetscInt off, node; 1905 for (v = 0; v < elem->numVerts; ++v) 1906 cone[v] = elem->nodes[v]; 1907 CHKERRQ(DMPlexReorderCell(cdm, cell, cone)); 1908 CHKERRQ(PetscSectionGetOffset(coordSection, cell, &off)); 1909 for (v = 0; v < elem->numVerts; ++v) 1910 for (node = cone[v], d = 0; d < coordDim; ++d) 1911 pointCoords[off++] = (PetscReal) coords[node*3+d]; 1912 } 1913 } 1914 } 1915 CHKERRQ(PetscMalloc1(numVerts, &nodeMap)); 1916 for (n = 0; n < numNodes; n++) 1917 if (mesh->vertexMap[n] >= 0) 1918 nodeMap[mesh->vertexMap[n]] = n; 1919 for (v = 0; v < numVerts; ++v) { 1920 PetscInt off, node = nodeMap[v]; 1921 CHKERRQ(PetscSectionGetOffset(coordSection, numCells + v, &off)); 1922 for (d = 0; d < coordDim; ++d) 1923 pointCoords[off+d] = (PetscReal) coords[node*3+d]; 1924 } 1925 CHKERRQ(PetscFree(nodeMap)); 1926 CHKERRQ(VecRestoreArray(coordinates, &pointCoords)); 1927 1928 } 1929 1930 CHKERRQ(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 1931 CHKERRQ(VecSetBlockSize(coordinates, coordDim)); 1932 CHKERRQ(DMSetCoordinatesLocal(*dm, coordinates)); 1933 CHKERRQ(VecDestroy(&coordinates)); 1934 CHKERRQ(DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL)); 1935 1936 CHKERRQ(GmshMeshDestroy(&mesh)); 1937 CHKERRQ(PetscBTDestroy(&periodicVerts)); 1938 CHKERRQ(PetscBTDestroy(&periodicCells)); 1939 1940 if (highOrder && project) { 1941 PetscFE fe; 1942 const char prefix[] = "dm_plex_gmsh_project_"; 1943 PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE; 1944 PetscDTNodeType nodeType = PETSCDTNODES_GAUSSJACOBI; 1945 1946 if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */ 1947 1948 CHKERRQ(GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe)); 1949 CHKERRQ(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view")); 1950 CHKERRQ(DMProjectCoordinates(*dm, fe)); 1951 CHKERRQ(PetscFEDestroy(&fe)); 1952 } 1953 1954 CHKERRQ(PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,NULL,NULL,NULL)); 1955 PetscFunctionReturn(0); 1956 } 1957