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