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