Actual source code: plexfluent.c
1: #define PETSC_DESIRE_FEATURE_TEST_MACROS /* for fileno() */
2: #define PETSCDM_DLL
3: #include <petsc/private/dmpleximpl.h>
5: /*@
6: DMPlexCreateFluentFromFile - Create a `DMPLEX` mesh from a Fluent mesh file
8: Collective
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
15: Output Parameter:
16: . dm - The `DM` object representing the mesh
18: Level: beginner
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;
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: }
37: static PetscErrorCode DMPlexCreateFluent_ReadString(PetscViewer viewer, char *buffer, char delim)
38: {
39: PetscInt ret, i = 0;
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: }
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;
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: }
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;
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));
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: }
107: static PetscErrorCode DMPlexCreateFluent_ReadSection(PetscViewer viewer, FluentSection *s)
108: {
109: char buffer[PETSC_MAX_PATH_LEN];
110: int snum;
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: }
123: if (s->index == 0) { /* Comment */
124: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
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");
131: } else if (s->index == 10 || s->index == 2010) { /* Vertices */
132: PetscInt numClosingParens = 0;
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;
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;
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];
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: }
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;
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: }
279: static PetscErrorCode ReorderPolygon(DM dm, PetscInt cell)
280: {
281: const PetscInt *cone, *ornt;
282: PetscInt coneSize, newCone[16], newOrnt[16];
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;
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;
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;
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: }
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];
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;
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;
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: }
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];
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;
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;
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;
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;
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]);
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: }
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: }
486: /*@C
487: DMPlexCreateFluent - Create a `DMPLEX` mesh from a Fluent mesh file <http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm>.
489: Collective
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
496: Output Parameter:
497: . dm - The `DM` object representing the mesh
499: Level: beginner
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;
527: PetscFunctionBegin;
528: PetscCallMPI(MPI_Comm_rank(comm, &rank));
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: }
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;
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;
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");
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));
648: if (rank == 0 && faces) {
649: PetscInt *cones;
651: PetscCall(DMPlexGetCones(*dm, &cones));
652: PetscCheck(numCellFaces < 16, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of cell faces %" PetscInt_FMT " exceeeds temporary storage", numCellFaces);
653: PetscCheck(numFaceVertices < 16, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of face vertices %" PetscInt_FMT " exceeeds 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];
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;
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: }
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];
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: }
722: { /* Create Face Sets label at all processes */
723: enum {
724: n = 1
725: };
726: PetscBool flag[n];
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: }
734: if (!interpolate) {
735: DM udm;
737: PetscCall(DMPlexUninterpolate(*dm, &udm));
738: PetscCall(DMDestroy(dm));
739: *dm = udm;
740: }
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));
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: }