1 #define PETSC_DESIRE_FEATURE_TEST_MACROS /* for fileno() */ 2 #define PETSCDM_DLL 3 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 4 5 /*@ 6 DMPlexCreateFluentFromFile - Create a `DMPLEX` mesh from a Fluent mesh file 7 8 Collective 9 10 Input Parameters: 11 + comm - The MPI communicator 12 . filename - Name of the Fluent mesh file 13 - interpolate - Create faces and edges in the mesh 14 15 Output Parameter: 16 . dm - The `DM` object representing the mesh 17 18 Level: beginner 19 20 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateFromFile()`, `DMPlexCreateFluent()`, `DMPlexCreate()` 21 @*/ 22 PetscErrorCode DMPlexCreateFluentFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 23 { 24 PetscViewer viewer; 25 26 PetscFunctionBegin; 27 /* Create file viewer and build plex */ 28 PetscCall(PetscViewerCreate(comm, &viewer)); 29 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII)); 30 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ)); 31 PetscCall(PetscViewerFileSetName(viewer, filename)); 32 PetscCall(DMPlexCreateFluent(comm, viewer, interpolate, dm)); 33 PetscCall(PetscViewerDestroy(&viewer)); 34 PetscFunctionReturn(PETSC_SUCCESS); 35 } 36 37 static PetscErrorCode DMPlexCreateFluent_ReadString(PetscViewer viewer, char *buffer, char delim) 38 { 39 PetscInt ret, i = 0; 40 41 PetscFunctionBegin; 42 do PetscCall(PetscViewerRead(viewer, &buffer[i++], 1, &ret, PETSC_CHAR)); 43 while (ret > 0 && buffer[i - 1] != '\0' && buffer[i - 1] != delim && i < PETSC_MAX_PATH_LEN - 1); 44 if (!ret) buffer[i - 1] = '\0'; 45 else buffer[i] = '\0'; 46 PetscCheck(i < PETSC_MAX_PATH_LEN - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Buffer overflow! This is not a valid Fluent file."); 47 PetscFunctionReturn(PETSC_SUCCESS); 48 } 49 50 static PetscErrorCode DMPlexCreateFluent_ReadValues(PetscViewer viewer, void *data, PetscInt count, PetscDataType dtype, PetscBool binary, PetscInt *numClosingParens) 51 { 52 int fdes = 0; 53 FILE *file; 54 PetscInt i; 55 56 PetscFunctionBegin; 57 *numClosingParens = 0; 58 if (binary) { 59 /* Extract raw file descriptor to read binary block */ 60 PetscCall(PetscViewerASCIIGetPointer(viewer, &file)); 61 PetscCall(PetscFFlush(file)); 62 fdes = fileno(file); 63 } 64 65 if (!binary && dtype == PETSC_INT) { 66 char cbuf[256]; 67 unsigned int ibuf; 68 int snum; 69 /* Parse hexadecimal ascii integers */ 70 for (i = 0; i < count; i++) { 71 size_t len; 72 73 PetscCall(PetscViewerRead(viewer, cbuf, 1, NULL, PETSC_STRING)); 74 snum = sscanf(cbuf, "%x", &ibuf); 75 PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 76 ((PetscInt *)data)[i] = (PetscInt)ibuf; 77 // Check for trailing parentheses 78 PetscCall(PetscStrlen(cbuf, &len)); 79 while (cbuf[len - 1] == ')' && len > 0) { 80 ++(*numClosingParens); 81 --len; 82 } 83 } 84 } else if (binary && dtype == PETSC_INT) { 85 /* Always read 32-bit ints and cast to PetscInt */ 86 int *ibuf; 87 PetscCall(PetscMalloc1(count, &ibuf)); 88 PetscCall(PetscBinaryRead(fdes, ibuf, count, NULL, PETSC_ENUM)); 89 PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, count)); 90 for (i = 0; i < count; i++) ((PetscInt *)data)[i] = ibuf[i]; 91 PetscCall(PetscFree(ibuf)); 92 93 } else if (binary && dtype == PETSC_SCALAR) { 94 float *fbuf; 95 /* Always read 32-bit floats and cast to PetscScalar */ 96 PetscCall(PetscMalloc1(count, &fbuf)); 97 PetscCall(PetscBinaryRead(fdes, fbuf, count, NULL, PETSC_FLOAT)); 98 PetscCall(PetscByteSwap(fbuf, PETSC_FLOAT, count)); 99 for (i = 0; i < count; i++) ((PetscScalar *)data)[i] = fbuf[i]; 100 PetscCall(PetscFree(fbuf)); 101 } else { 102 PetscCall(PetscViewerASCIIRead(viewer, data, count, NULL, dtype)); 103 } 104 PetscFunctionReturn(PETSC_SUCCESS); 105 } 106 107 static PetscErrorCode DMPlexCreateFluent_ReadSection(PetscViewer viewer, FluentSection *s) 108 { 109 char buffer[PETSC_MAX_PATH_LEN]; 110 int snum; 111 112 PetscFunctionBegin; 113 /* Fast-forward to next section and derive its index */ 114 PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '(')); 115 PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ' ')); 116 snum = sscanf(buffer, "%d", &s->index); 117 /* If we can't match an index return -1 to signal end-of-file */ 118 if (snum < 1) { 119 s->index = -1; 120 PetscFunctionReturn(PETSC_SUCCESS); 121 } 122 123 if (s->index == 0) { /* Comment */ 124 PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')')); 125 126 } else if (s->index == 2) { /* Dimension */ 127 PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')')); 128 snum = sscanf(buffer, "%d", &s->nd); 129 PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 130 131 } else if (s->index == 10 || s->index == 2010) { /* Vertices */ 132 PetscInt numClosingParens = 0; 133 134 PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')')); 135 snum = sscanf(buffer, "(%x %x %x %d %d)", &s->zoneID, &s->first, &s->last, &s->type, &s->nd); 136 PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 137 if (s->zoneID > 0) { 138 PetscInt numCoords = s->last - s->first + 1; 139 PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '(')); 140 PetscCall(PetscMalloc1(s->nd * numCoords, (PetscScalar **)&s->data)); 141 PetscCall(DMPlexCreateFluent_ReadValues(viewer, s->data, s->nd * numCoords, PETSC_SCALAR, s->index == 2010 ? PETSC_TRUE : PETSC_FALSE, &numClosingParens)); 142 if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')')); 143 else --numClosingParens; 144 } 145 if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')')); 146 else --numClosingParens; 147 PetscCheck(!numClosingParens, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 148 } else if (s->index == 12 || s->index == 2012) { /* Cells */ 149 PetscInt numClosingParens = 0; 150 151 PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')')); 152 snum = sscanf(buffer, "(%x", &s->zoneID); 153 PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 154 if (s->zoneID == 0) { /* Header section */ 155 snum = sscanf(buffer, "(%x %x %x %d)", &s->zoneID, &s->first, &s->last, &s->nd); 156 PetscCheck(snum == 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 157 } else { /* Data section */ 158 snum = sscanf(buffer, "(%x %x %x %d %d)", &s->zoneID, &s->first, &s->last, &s->type, &s->nd); 159 PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 160 if (s->nd == 0) { 161 /* Read cell type definitions for mixed cells */ 162 PetscInt numCells = s->last - s->first + 1; 163 PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '(')); 164 PetscCall(PetscMalloc1(numCells, (PetscInt **)&s->data)); 165 PetscCall(DMPlexCreateFluent_ReadValues(viewer, s->data, numCells, PETSC_INT, s->index == 2012 ? PETSC_TRUE : PETSC_FALSE, &numClosingParens)); 166 PetscCall(PetscFree(s->data)); 167 if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')')); 168 else --numClosingParens; 169 } 170 } 171 if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')')); 172 else --numClosingParens; 173 PetscCheck(!numClosingParens, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 174 } else if (s->index == 13 || s->index == 2013) { /* Faces */ 175 PetscInt numClosingParens = 0; 176 177 PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')')); 178 snum = sscanf(buffer, "(%x", &s->zoneID); 179 PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 180 if (s->zoneID == 0) { /* Header section */ 181 snum = sscanf(buffer, "(%x %x %x %d)", &s->zoneID, &s->first, &s->last, &s->nd); 182 PetscCheck(snum == 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 183 } else { /* Data section */ 184 PetscInt f, numEntries, numFaces; 185 snum = sscanf(buffer, "(%x %x %x %d %d)", &s->zoneID, &s->first, &s->last, &s->type, &s->nd); 186 PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 187 PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '(')); 188 switch (s->nd) { 189 case 0: 190 numEntries = PETSC_DETERMINE; 191 break; 192 case 2: 193 numEntries = 2 + 2; 194 break; /* linear */ 195 case 3: 196 numEntries = 2 + 3; 197 break; /* triangular */ 198 case 4: 199 numEntries = 2 + 4; 200 break; /* quadrilateral */ 201 default: 202 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file"); 203 } 204 numFaces = s->last - s->first + 1; 205 if (numEntries != PETSC_DETERMINE) { 206 /* Allocate space only if we already know the size of the block */ 207 PetscCall(PetscMalloc1(numEntries * numFaces, (PetscInt **)&s->data)); 208 } 209 for (f = 0; f < numFaces; f++) { 210 if (s->nd == 0) { 211 /* Determine the size of the block for "mixed" facets */ 212 PetscInt numFaceVert = 0; 213 PetscCall(DMPlexCreateFluent_ReadValues(viewer, &numFaceVert, 1, PETSC_INT, s->index == 2013 ? PETSC_TRUE : PETSC_FALSE, &numClosingParens)); 214 if (numEntries == PETSC_DETERMINE) { 215 numEntries = numFaceVert + 2; 216 PetscCall(PetscMalloc1(numEntries * numFaces, (PetscInt **)&s->data)); 217 } else { 218 PetscCheck(numEntries == numFaceVert + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for mixed faces in Fluent files"); 219 } 220 } 221 PetscCall(DMPlexCreateFluent_ReadValues(viewer, &(((PetscInt *)s->data)[f * numEntries]), numEntries, PETSC_INT, s->index == 2013 ? PETSC_TRUE : PETSC_FALSE, &numClosingParens)); 222 } 223 PetscCall(PetscMPIIntCast(numEntries - 2, &s->nd)); 224 if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')')); 225 else --numClosingParens; 226 } 227 if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')')); 228 else --numClosingParens; 229 PetscCheck(!numClosingParens, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 230 } else if (s->index == 39) { /* Label information */ 231 char labelName[PETSC_MAX_PATH_LEN]; 232 char caseName[PETSC_MAX_PATH_LEN]; 233 234 PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')')); 235 snum = sscanf(buffer, "(%u %s %s %d)", &s->zoneID, caseName, labelName, &s->nd); 236 PetscCheck(snum == 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file: %d", snum); 237 PetscInt depth = 1; 238 do { 239 /* Match parentheses when parsing unknown sections */ 240 do PetscCall(PetscViewerRead(viewer, &buffer[0], 1, NULL, PETSC_CHAR)); 241 while (buffer[0] != '(' && buffer[0] != ')'); 242 if (buffer[0] == '(') depth++; 243 if (buffer[0] == ')') depth--; 244 } while (depth > 0); 245 PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '\n')); 246 PetscCall(PetscStrallocpy(labelName, (char **)&s->data)); 247 PetscCall(PetscInfo((PetscObject)viewer, "CASE: Zone ID %u is label %s\n", s->zoneID, labelName)); 248 } else { /* Unknown section type */ 249 PetscInt depth = 1; 250 do { 251 /* Match parentheses when parsing unknown sections */ 252 do PetscCall(PetscViewerRead(viewer, &buffer[0], 1, NULL, PETSC_CHAR)); 253 while (buffer[0] != '(' && buffer[0] != ')'); 254 if (buffer[0] == '(') depth++; 255 if (buffer[0] == ')') depth--; 256 } while (depth > 0); 257 PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '\n')); 258 } 259 PetscFunctionReturn(PETSC_SUCCESS); 260 } 261 262 // Inserts point `face` with orientation `ornt` into the cone of point `cell` at position `c`, which is the first empty slot 263 static PetscErrorCode InsertFace(DM dm, PetscInt cell, PetscInt face, PetscInt ornt) 264 { 265 const PetscInt *cone; 266 PetscInt coneSize, c; 267 268 PetscFunctionBegin; 269 PetscCall(DMPlexGetCone(dm, cell, &cone)); 270 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 271 for (c = 0; c < coneSize; ++c) 272 if (cone[c] < 0) break; 273 PetscCheck(c < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be inserted in cone of cell %" PetscInt_FMT, face, cell); 274 PetscCall(DMPlexInsertCone(dm, cell, c, face)); 275 PetscCall(DMPlexInsertConeOrientation(dm, cell, c, ornt)); 276 PetscFunctionReturn(PETSC_SUCCESS); 277 } 278 279 static PetscErrorCode ReorderPolygon(DM dm, PetscInt cell) 280 { 281 const PetscInt *cone, *ornt; 282 PetscInt coneSize, newCone[16], newOrnt[16]; 283 284 PetscFunctionBegin; 285 PetscCall(DMPlexGetOrientedCone(dm, cell, &cone, &ornt)); 286 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 287 newCone[0] = cone[0]; 288 newOrnt[0] = ornt[0]; 289 for (PetscInt c = 1; c < coneSize; ++c) { 290 const PetscInt *fcone; 291 PetscInt firstVertex, lastVertex, c2; 292 293 PetscCall(DMPlexGetCone(dm, newCone[c - 1], &fcone)); 294 lastVertex = newOrnt[c - 1] ? fcone[0] : fcone[1]; 295 for (c2 = 0; c2 < coneSize; ++c2) { 296 const PetscInt *fcone2; 297 298 PetscCall(DMPlexGetCone(dm, cone[c2], &fcone2)); 299 firstVertex = ornt[c2] ? fcone2[1] : fcone2[0]; 300 if (lastVertex == firstVertex) { 301 // Point `cell` matched point `lastVertex` on face `cone[c2]` with orientation `ornt[c2]` 302 break; 303 } 304 } 305 PetscCheck(c2 < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " could not find a face match as position %" PetscInt_FMT, cell, c); 306 newCone[c] = cone[c2]; 307 newOrnt[c] = ornt[c2]; 308 } 309 { 310 const PetscInt *fcone, *fcone2; 311 PetscInt vertex, vertex2; 312 313 PetscCall(DMPlexGetCone(dm, newCone[coneSize - 1], &fcone)); 314 PetscCall(DMPlexGetCone(dm, newCone[0], &fcone2)); 315 vertex = newOrnt[coneSize - 1] ? fcone[0] : fcone[1]; 316 vertex2 = newOrnt[0] ? fcone2[1] : fcone2[0]; 317 PetscCheck(vertex == vertex2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " did not match at the endpoint", cell); 318 } 319 PetscCall(DMPlexSetCone(dm, cell, newCone)); 320 PetscCall(DMPlexSetConeOrientation(dm, cell, newOrnt)); 321 PetscCall(DMPlexRestoreOrientedCone(dm, cell, &cone, &ornt)); 322 PetscFunctionReturn(PETSC_SUCCESS); 323 } 324 325 static PetscErrorCode ReorderTetrahedron(DM dm, PetscInt cell) 326 { 327 const PetscInt *cone, *ornt, *fcone, *fornt, *farr, faces[4] = {0, 1, 3, 2}; 328 PetscInt newCone[16], newOrnt[16]; 329 330 PetscFunctionBegin; 331 PetscCall(DMPlexGetOrientedCone(dm, cell, &cone, &ornt)); 332 newCone[0] = cone[0]; 333 newOrnt[0] = ornt[0]; 334 PetscCall(DMPlexGetOrientedCone(dm, newCone[0], &fcone, &fornt)); 335 farr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, newOrnt[0]); 336 // Loop over each edge in the initial triangle 337 for (PetscInt e = 0; e < 3; ++e) { 338 const PetscInt edge = fcone[farr[e * 2 + 0]], eornt = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[e * 2 + 1], fornt[farr[e * 2 + 0]]); 339 PetscInt c; 340 341 // Loop over each remaining face in the tetrahedron 342 // On face `newCone[0]`, trying to match edge `edge` with final orientation `eornt` to an edge on another face 343 for (c = 1; c < 4; ++c) { 344 const PetscInt *fcone2, *fornt2, *farr2; 345 PetscInt c2; 346 347 // Checking face `cone[c]` with orientation `ornt[c]` 348 PetscCall(DMPlexGetOrientedCone(dm, cone[c], &fcone2, &fornt2)); 349 farr2 = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, ornt[c]); 350 // Check for edge 351 for (c2 = 0; c2 < 3; ++c2) { 352 const PetscInt edge2 = fcone2[farr2[c2 * 2 + 0]], eornt2 = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[c2 * 2 + 1], fornt2[farr2[c2 * 2 + 0]]); 353 // Trying to match edge `edge2` with final orientation `eornt2` 354 if (edge == edge2) { 355 PetscCheck(eornt == -(eornt2 + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " found twice with the same orientation", edge); 356 // Matched face `newCone[0]` with orientation `newOrnt[0]` to face `cone[c]` with orientation `ornt[c]` along edge `edge` 357 break; 358 } 359 } 360 if (c2 < 3) { 361 newCone[faces[e + 1]] = cone[c]; 362 // Compute new orientation of face based on which edge was matched (only the first edge matches a side different from 0) 363 newOrnt[faces[e + 1]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_TRIANGLE, !e ? (c2 + 1) % 3 : c2, ornt[c]); 364 break; 365 } 366 } 367 PetscCheck(c < 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " could not find a face match for edge %" PetscInt_FMT, cell, e); 368 } 369 PetscCall(DMPlexRestoreOrientedCone(dm, newCone[0], &fcone, &fornt)); 370 PetscCall(DMPlexSetCone(dm, cell, newCone)); 371 PetscCall(DMPlexSetConeOrientation(dm, cell, newOrnt)); 372 PetscCall(DMPlexRestoreOrientedCone(dm, cell, &cone, &ornt)); 373 PetscFunctionReturn(PETSC_SUCCESS); 374 } 375 376 static PetscErrorCode ReorderHexahedron(DM dm, PetscInt cell) 377 { 378 const PetscInt *cone, *ornt, *fcone, *fornt, *farr; 379 const PetscInt faces[6] = {0, 5, 3, 4, 2, 1}; 380 PetscInt used[6] = {1, 0, 0, 0, 0, 0}; 381 PetscInt newCone[16], newOrnt[16]; 382 383 PetscFunctionBegin; 384 PetscCall(DMPlexGetOrientedCone(dm, cell, &cone, &ornt)); 385 newCone[0] = cone[0]; 386 newOrnt[0] = ornt[0]; 387 PetscCall(DMPlexGetOrientedCone(dm, newCone[0], &fcone, &fornt)); 388 farr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, newOrnt[0]); 389 // Loop over each edge in the initial quadrilateral 390 for (PetscInt e = 0; e < 4; ++e) { 391 const PetscInt edge = fcone[farr[e * 2 + 0]], eornt = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[e * 2 + 1], fornt[farr[e * 2 + 0]]); 392 PetscInt c; 393 394 // Loop over each remaining face in the hexahedron 395 // On face `newCone[0]`, trying to match edge `edge` with final orientation `eornt` to an edge on another face 396 for (c = 1; c < 6; ++c) { 397 const PetscInt *fcone2, *fornt2, *farr2; 398 PetscInt c2; 399 400 // Checking face `cone[c]` with orientation `ornt[c]` 401 PetscCall(DMPlexGetOrientedCone(dm, cone[c], &fcone2, &fornt2)); 402 farr2 = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, ornt[c]); 403 // Check for edge 404 for (c2 = 0; c2 < 4; ++c2) { 405 const PetscInt edge2 = fcone2[farr2[c2 * 2 + 0]], eornt2 = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr2[c2 * 2 + 1], fornt2[farr2[c2 * 2 + 0]]); 406 // Trying to match edge `edge2` with final orientation `eornt2` 407 if (edge == edge2) { 408 PetscCheck(eornt == -(eornt2 + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " found twice with the same orientation", edge); 409 // Matched face `newCone[0]` with orientation `newOrnt[0]` to face `cone[c]` with orientation `ornt[c]` along edge `edge` 410 break; 411 } 412 } 413 if (c2 < 4) { 414 used[c] = 1; 415 newCone[faces[e + 1]] = cone[c]; 416 // Compute new orientation of face based on which edge was matched (only the first edge matches a side different from 0) 417 newOrnt[faces[e + 1]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_QUADRILATERAL, !e ? (c2 + 1) % 4 : c2, ornt[c]); 418 break; 419 } 420 } 421 PetscCheck(c < 6, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " could not find a face match for edge %" PetscInt_FMT, cell, e); 422 } 423 PetscCall(DMPlexRestoreOrientedCone(dm, newCone[0], &fcone, &fornt)); 424 // Add last face 425 { 426 PetscInt c, c2; 427 428 for (c = 1; c < 6; ++c) 429 if (!used[c]) break; 430 PetscCheck(c < 6, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " could not find an available face", cell); 431 // Match first edge to 3rd edge in newCone[2] 432 { 433 const PetscInt *fcone2, *fornt2, *farr2; 434 435 PetscCall(DMPlexGetOrientedCone(dm, newCone[2], &fcone, &fornt)); 436 farr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, newOrnt[2]); 437 PetscCall(DMPlexGetOrientedCone(dm, cone[c], &fcone2, &fornt2)); 438 farr2 = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, ornt[c]); 439 440 const PetscInt e = 2; 441 const PetscInt edge = fcone[farr[e * 2 + 0]], eornt = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[e * 2 + 1], fornt[farr[e * 2 + 0]]); 442 // Trying to match edge `edge` with final orientation `eornt` of face `newCone[2]` to some edge of face `cone[c]` with orientation `ornt[c]` 443 for (c2 = 0; c2 < 4; ++c2) { 444 const PetscInt edge2 = fcone2[farr2[c2 * 2 + 0]], eornt2 = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr2[c2 * 2 + 1], fornt2[farr2[c2 * 2 + 0]]); 445 // Trying to match edge `edge2` with final orientation `eornt2` 446 if (edge == edge2) { 447 PetscCheck(eornt == -(eornt2 + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " found twice with the same orientation", edge); 448 // Matched face `newCone[2]` with orientation `newOrnt[2]` to face `cone[c]` with orientation `ornt[c]` along edge `edge` 449 break; 450 } 451 } 452 PetscCheck(c2 < 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not fit last face in"); 453 } 454 newCone[faces[5]] = cone[c]; 455 // Compute new orientation of face based on which edge was matched 456 newOrnt[faces[5]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_QUADRILATERAL, c2, ornt[c]); 457 PetscCall(DMPlexRestoreOrientedCone(dm, newCone[0], &fcone, &fornt)); 458 } 459 PetscCall(DMPlexSetCone(dm, cell, newCone)); 460 PetscCall(DMPlexSetConeOrientation(dm, cell, newOrnt)); 461 PetscCall(DMPlexRestoreOrientedCone(dm, cell, &cone, &ornt)); 462 PetscFunctionReturn(PETSC_SUCCESS); 463 } 464 465 static PetscErrorCode ReorderCell(DM dm, PetscInt cell, DMPolytopeType ct) 466 { 467 PetscFunctionBegin; 468 switch (ct) { 469 case DM_POLYTOPE_TRIANGLE: 470 case DM_POLYTOPE_QUADRILATERAL: 471 PetscCall(ReorderPolygon(dm, cell)); 472 break; 473 case DM_POLYTOPE_TETRAHEDRON: 474 PetscCall(ReorderTetrahedron(dm, cell)); 475 break; 476 case DM_POLYTOPE_HEXAHEDRON: 477 PetscCall(ReorderHexahedron(dm, cell)); 478 break; 479 default: 480 PetscCheck(0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Celltype %s is unsupported", DMPolytopeTypes[ct]); 481 break; 482 } 483 PetscFunctionReturn(PETSC_SUCCESS); 484 } 485 486 /*@C 487 DMPlexCreateFluent - Create a `DMPLEX` mesh from a Fluent mesh file <http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm>. 488 489 Collective 490 491 Input Parameters: 492 + comm - The MPI communicator 493 . viewer - The `PetscViewer` associated with a Fluent mesh file 494 - interpolate - Create faces and edges in the mesh 495 496 Output Parameter: 497 . dm - The `DM` object representing the mesh 498 499 Level: beginner 500 501 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()` 502 @*/ 503 PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm) 504 { 505 PetscInt dim = PETSC_DETERMINE; 506 PetscInt numCells = 0; 507 PetscInt numVertices = 0; 508 PetscInt numCellFaces = PETSC_DETERMINE; 509 DMPolytopeType ct = DM_POLYTOPE_UNKNOWN_CELL; 510 PetscInt numFaces = 0; 511 PetscInt numFaceEntries = PETSC_DETERMINE; 512 PetscInt numFaceVertices = PETSC_DETERMINE; 513 PetscInt *faces = NULL; 514 PetscInt *cellVertices = NULL; 515 unsigned int *faceZoneIDs = NULL; 516 DMLabel faceSets = NULL; 517 DMLabel *zoneLabels = NULL; 518 const char **zoneNames = NULL; 519 unsigned int maxZoneID = 0; 520 PetscScalar *coordsIn = NULL; 521 PetscScalar *coords; 522 PetscSection coordSection; 523 Vec coordinates; 524 PetscInt coordSize, f; 525 PetscMPIInt rank; 526 527 PetscFunctionBegin; 528 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 529 530 if (rank == 0) { 531 FluentSection s; 532 numFaces = PETSC_DETERMINE; 533 do { 534 PetscCall(DMPlexCreateFluent_ReadSection(viewer, &s)); 535 if (s.index == 2) { /* Dimension */ 536 dim = s.nd; 537 PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found dimension: %" PetscInt_FMT "\n", dim)); 538 } else if (s.index == 10 || s.index == 2010) { /* Vertices */ 539 if (s.zoneID == 0) { 540 numVertices = s.last; 541 PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found number of vertices: %" PetscInt_FMT "\n", numVertices)); 542 } else { 543 PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found vertex coordinates\n")); 544 PetscCheck(!coordsIn, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files"); 545 coordsIn = (PetscScalar *)s.data; 546 } 547 548 } else if (s.index == 12 || s.index == 2012) { /* Cells */ 549 if (s.zoneID == 0) { 550 numCells = s.last; 551 PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found number of cells %" PetscInt_FMT "\n", numCells)); 552 } else { 553 switch (s.nd) { 554 case 0: 555 numCellFaces = PETSC_DETERMINE; 556 ct = DM_POLYTOPE_POINT; 557 break; 558 case 1: 559 numCellFaces = 3; 560 ct = DM_POLYTOPE_TRIANGLE; 561 break; 562 case 2: 563 numCellFaces = 4; 564 ct = DM_POLYTOPE_TETRAHEDRON; 565 break; 566 case 3: 567 numCellFaces = 4; 568 ct = DM_POLYTOPE_QUADRILATERAL; 569 break; 570 case 4: 571 numCellFaces = 6; 572 ct = DM_POLYTOPE_HEXAHEDRON; 573 break; 574 case 5: 575 numCellFaces = 5; 576 ct = DM_POLYTOPE_PYRAMID; 577 break; 578 case 6: 579 numCellFaces = 5; 580 ct = DM_POLYTOPE_TRI_PRISM; 581 break; 582 default: 583 numCellFaces = PETSC_DETERMINE; 584 } 585 PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found number of cell faces %" PetscInt_FMT "\n", numCellFaces)); 586 } 587 } else if (s.index == 13 || s.index == 2013) { /* Facets */ 588 if (s.zoneID == 0) { /* Header section */ 589 numFaces = (PetscInt)(s.last - s.first + 1); 590 if (s.nd == 0 || s.nd == 5) numFaceVertices = PETSC_DETERMINE; 591 else numFaceVertices = s.nd; 592 PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found number of faces %" PetscInt_FMT " face vertices: %" PetscInt_FMT "\n", numFaces, numFaceVertices)); 593 } else { /* Data section */ 594 unsigned int z; 595 596 PetscCheck(numFaceVertices == PETSC_DETERMINE || s.nd == numFaceVertices, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed facets in Fluent files are not supported"); 597 PetscCheck(numFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file"); 598 if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd; 599 numFaceEntries = numFaceVertices + 2; 600 if (!faces) PetscCall(PetscMalloc1(numFaces * numFaceEntries, &faces)); 601 if (!faceZoneIDs) PetscCall(PetscMalloc1(numFaces, &faceZoneIDs)); 602 PetscCall(PetscMemcpy(&faces[(s.first - 1) * numFaceEntries], s.data, (s.last - s.first + 1) * numFaceEntries * sizeof(PetscInt))); 603 /* Record the zoneID for each face set */ 604 for (z = s.first - 1; z < s.last; z++) faceZoneIDs[z] = s.zoneID; 605 PetscCall(PetscFree(s.data)); 606 } 607 } else if (s.index == 39) { /* Label information */ 608 if (s.zoneID >= maxZoneID) { 609 DMLabel *tmpL; 610 const char **tmp; 611 unsigned int newmax = maxZoneID + 1; 612 613 while (newmax < s.zoneID + 1) newmax *= 2; 614 PetscCall(PetscCalloc2(newmax, &tmp, newmax, &tmpL)); 615 for (PetscInt i = 0; i < (PetscInt)maxZoneID; ++i) { 616 tmp[i] = zoneNames[i]; 617 tmpL[i] = zoneLabels[i]; 618 } 619 maxZoneID = newmax; 620 PetscCall(PetscFree2(zoneNames, zoneLabels)); 621 zoneNames = tmp; 622 zoneLabels = tmpL; 623 } 624 zoneNames[s.zoneID] = (const char *)s.data; 625 } 626 } while (s.index >= 0); 627 } 628 PetscCallMPI(MPI_Bcast(&dim, 1, MPIU_INT, 0, comm)); 629 PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension"); 630 631 /* Allocate cell-vertex mesh */ 632 PetscCall(DMCreate(comm, dm)); 633 PetscCall(DMSetType(*dm, DMPLEX)); 634 PetscCall(DMSetDimension(*dm, dim)); 635 PetscCall(DMPlexSetChart(*dm, 0, numCells + numFaces + numVertices)); 636 if (rank == 0) { 637 PetscCheck(numCells >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown number of cells in Fluent file"); 638 /* If no cell type was given we assume simplices */ 639 if (numCellFaces == PETSC_DETERMINE) { 640 numCellFaces = numFaceVertices + 1; 641 ct = numCellFaces == 3 ? DM_POLYTOPE_TRIANGLE : (numCellFaces == 4 ? DM_POLYTOPE_TETRAHEDRON : DM_POLYTOPE_UNKNOWN_CELL); 642 } 643 for (PetscInt c = 0; c < numCells; ++c) PetscCall(DMPlexSetConeSize(*dm, c, numCellFaces)); 644 for (PetscInt f = 0; f < numFaces; ++f) PetscCall(DMPlexSetConeSize(*dm, f + numCells + numVertices, numFaceVertices)); 645 } 646 PetscCall(DMSetUp(*dm)); 647 648 if (rank == 0 && faces) { 649 PetscInt *cones; 650 651 PetscCall(DMPlexGetCones(*dm, &cones)); 652 PetscCheck(numCellFaces < 16, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of cell faces %" PetscInt_FMT " exceeds temporary storage", numCellFaces); 653 PetscCheck(numFaceVertices < 16, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of face vertices %" PetscInt_FMT " exceeds temporary storage", numFaceVertices); 654 for (PetscInt c = 0; c < numCells * numCellFaces; ++c) cones[c] = -1; 655 for (PetscInt f = 0; f < numFaces; f++) { 656 const PetscInt cl = faces[f * numFaceEntries + numFaceVertices] - 1; 657 const PetscInt cr = faces[f * numFaceEntries + numFaceVertices + 1] - 1; 658 const PetscInt face = f + numCells + numVertices; 659 const PetscInt *fc = &faces[f * numFaceEntries]; 660 PetscInt fcone[16]; 661 662 // How could Fluent define the outward normal differently? Is there no end to the pain? 663 if (dim == 3) { 664 if (cl >= 0) PetscCall(InsertFace(*dm, cl, face, -1)); 665 if (cr >= 0) PetscCall(InsertFace(*dm, cr, face, 0)); 666 } else { 667 if (cl >= 0) PetscCall(InsertFace(*dm, cl, face, 0)); 668 if (cr >= 0) PetscCall(InsertFace(*dm, cr, face, -1)); 669 } 670 for (PetscInt v = 0; v < numFaceVertices; ++v) fcone[v] = fc[v] + numCells - 1; 671 PetscCall(DMPlexSetCone(*dm, face, fcone)); 672 } 673 } 674 PetscCall(DMPlexSymmetrize(*dm)); 675 PetscCall(DMPlexStratify(*dm)); 676 if (dim == 3) { 677 DM idm; 678 679 PetscCall(DMCreate(PetscObjectComm((PetscObject)*dm), &idm)); 680 PetscCall(DMSetType(idm, DMPLEX)); 681 PetscCall(DMSetDimension(idm, dim)); 682 PetscCall(DMPlexInterpolateFaces_Internal(*dm, 1, idm)); 683 PetscCall(DMDestroy(dm)); 684 *dm = idm; 685 } 686 PetscCall(DMViewFromOptions(*dm, NULL, "-cas_dm_view")); 687 if (rank == 0 && faces) { 688 for (PetscInt c = 0; c < numCells; ++c) PetscCall(ReorderCell(*dm, c, ct)); 689 } 690 691 if (rank == 0 && faces) { 692 PetscInt fi, joinSize, meetSize, *fverts, cells[2]; 693 const PetscInt *join, *meet; 694 PetscCall(PetscMalloc1(numFaceVertices, &fverts)); 695 /* Mark facets by finding the full join of all adjacent vertices */ 696 for (f = 0; f < numFaces; f++) { 697 const PetscInt cl = faces[f * numFaceEntries + numFaceVertices] - 1; 698 const PetscInt cr = faces[f * numFaceEntries + numFaceVertices + 1] - 1; 699 const PetscInt id = (PetscInt)faceZoneIDs[f]; 700 701 if (cl > 0 && cr > 0) { 702 /* If we know both adjoining cells we can use a single-level meet */ 703 cells[0] = cl; 704 cells[1] = cr; 705 PetscCall(DMPlexGetMeet(*dm, 2, cells, &meetSize, &meet)); 706 PetscCheck(meetSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %" PetscInt_FMT " cells: %" PetscInt_FMT ", %" PetscInt_FMT, f, cl, cr); 707 PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", meet[0], id)); 708 if (zoneNames && zoneNames[id]) PetscCall(DMSetLabelValue_Fast(*dm, &zoneLabels[id], zoneNames[id], meet[0], 1)); 709 PetscCall(DMPlexRestoreMeet(*dm, numFaceVertices, fverts, &meetSize, &meet)); 710 } else { 711 for (fi = 0; fi < numFaceVertices; fi++) fverts[fi] = faces[f * numFaceEntries + fi] + numCells - 1; 712 PetscCall(DMPlexGetFullJoin(*dm, numFaceVertices, fverts, &joinSize, &join)); 713 PetscCheck(joinSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %" PetscInt_FMT, f); 714 PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], id)); 715 if (zoneNames && zoneNames[id]) PetscCall(DMSetLabelValue_Fast(*dm, &zoneLabels[id], zoneNames[id], join[0], 1)); 716 PetscCall(DMPlexRestoreJoin(*dm, numFaceVertices, fverts, &joinSize, &join)); 717 } 718 } 719 PetscCall(PetscFree(fverts)); 720 } 721 722 { /* Create Face Sets label at all processes */ 723 enum { 724 n = 1 725 }; 726 PetscBool flag[n]; 727 728 flag[0] = faceSets ? PETSC_TRUE : PETSC_FALSE; 729 PetscCallMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm)); 730 if (flag[0]) PetscCall(DMCreateLabel(*dm, "Face Sets")); 731 // TODO Code to create all the zone labels on each process 732 } 733 734 if (!interpolate) { 735 DM udm; 736 737 PetscCall(DMPlexUninterpolate(*dm, &udm)); 738 PetscCall(DMDestroy(dm)); 739 *dm = udm; 740 } 741 742 /* Read coordinates */ 743 PetscCall(DMGetCoordinateSection(*dm, &coordSection)); 744 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 745 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dim)); 746 PetscCall(PetscSectionSetChart(coordSection, numCells, numCells + numVertices)); 747 for (PetscInt v = numCells; v < numCells + numVertices; ++v) { 748 PetscCall(PetscSectionSetDof(coordSection, v, dim)); 749 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dim)); 750 } 751 PetscCall(PetscSectionSetUp(coordSection)); 752 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 753 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 754 PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates")); 755 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 756 PetscCall(VecSetType(coordinates, VECSTANDARD)); 757 PetscCall(VecGetArray(coordinates, &coords)); 758 if (rank == 0 && coordsIn) { 759 for (PetscInt v = 0; v < numVertices; ++v) { 760 for (PetscInt d = 0; d < dim; ++d) coords[v * dim + d] = coordsIn[v * dim + d]; 761 } 762 } 763 PetscCall(VecRestoreArray(coordinates, &coords)); 764 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 765 PetscCall(VecDestroy(&coordinates)); 766 767 if (rank == 0) { 768 PetscCall(PetscFree(cellVertices)); 769 PetscCall(PetscFree(faces)); 770 PetscCall(PetscFree(faceZoneIDs)); 771 PetscCall(PetscFree(coordsIn)); 772 if (zoneNames) 773 for (PetscInt i = 0; i < (PetscInt)maxZoneID; ++i) PetscCall(PetscFree(zoneNames[i])); 774 PetscCall(PetscFree2(zoneNames, zoneLabels)); 775 } 776 PetscFunctionReturn(PETSC_SUCCESS); 777 } 778