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