Actual source code: plexfluent.c
1: #define PETSC_DESIRE_FEATURE_TEST_MACROS /* for fileno() */
2: #include <petsc/private/dmpleximpl.h>
4: /* Utility struct to store the contents of a Fluent file in memory */
5: typedef struct {
6: int index; /* Type of section */
7: unsigned int zoneID;
8: unsigned int first;
9: unsigned int last;
10: int type;
11: int nd; /* Either ND or element-type */
12: void *data;
13: } FluentSection;
15: /*@
16: DMPlexCreateFluentFromFile - Create a `DMPLEX` mesh from a Fluent mesh file
18: Collective
20: Input Parameters:
21: + comm - The MPI communicator
22: . filename - Name of the Fluent mesh file
23: - interpolate - Create faces and edges in the mesh
25: Output Parameter:
26: . dm - The `DM` object representing the mesh
28: Level: beginner
30: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateFromFile()`, `DMPlexCreateFluent()`, `DMPlexCreate()`
31: @*/
32: PetscErrorCode DMPlexCreateFluentFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
33: {
34: PetscViewer viewer;
36: PetscFunctionBegin;
37: /* Create file viewer and build plex */
38: PetscCall(PetscViewerCreate(comm, &viewer));
39: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII));
40: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ));
41: PetscCall(PetscViewerFileSetName(viewer, filename));
42: PetscCall(DMPlexCreateFluent(comm, viewer, interpolate, dm));
43: PetscCall(PetscViewerDestroy(&viewer));
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: static PetscErrorCode DMPlexCreateFluent_ReadString(PetscViewer viewer, char *buffer, char delim)
48: {
49: PetscInt ret, i = 0;
51: PetscFunctionBegin;
52: do PetscCall(PetscViewerRead(viewer, &buffer[i++], 1, &ret, PETSC_CHAR));
53: while (ret > 0 && buffer[i - 1] != '\0' && buffer[i - 1] != delim && i < PETSC_MAX_PATH_LEN - 1);
54: if (!ret) buffer[i - 1] = '\0';
55: else buffer[i] = '\0';
56: PetscCheck(i < PETSC_MAX_PATH_LEN - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Buffer overflow! This is not a valid Fluent file.");
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: static PetscErrorCode DMPlexCreateFluent_ReadValues(PetscViewer viewer, void *data, PetscInt count, PetscDataType dtype, PetscBool binary, PetscInt *numClosingParens)
61: {
62: int fdes = 0;
63: FILE *file;
65: PetscFunctionBegin;
66: *numClosingParens = 0;
67: if (binary) {
68: /* Extract raw file descriptor to read binary block */
69: PetscCall(PetscViewerASCIIGetPointer(viewer, &file));
70: PetscCall(PetscFFlush(file));
71: fdes = fileno(file);
72: }
74: if (!binary && dtype == PETSC_INT) {
75: char cbuf[256];
76: unsigned int ibuf;
77: int snum;
78: /* Parse hexadecimal ascii integers */
79: for (PetscInt i = 0; i < count; i++) {
80: size_t len;
82: PetscCall(PetscViewerRead(viewer, cbuf, 1, NULL, PETSC_STRING));
83: snum = sscanf(cbuf, "%x", &ibuf);
84: PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
85: ((PetscInt *)data)[i] = (PetscInt)ibuf;
86: // Check for trailing parentheses
87: PetscCall(PetscStrlen(cbuf, &len));
88: while (cbuf[len - 1] == ')' && len > 0) {
89: ++(*numClosingParens);
90: --len;
91: }
92: }
93: } else if (binary && dtype == PETSC_INT) {
94: /* Always read 32-bit ints and cast to PetscInt */
95: int *ibuf;
96: PetscCall(PetscMalloc1(count, &ibuf));
97: PetscCall(PetscBinaryRead(fdes, ibuf, count, NULL, PETSC_ENUM));
98: PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, count));
99: for (PetscInt i = 0; i < count; i++) ((PetscInt *)data)[i] = ibuf[i];
100: PetscCall(PetscFree(ibuf));
102: } else if (binary && dtype == PETSC_SCALAR) {
103: float *fbuf;
104: /* Always read 32-bit floats and cast to PetscScalar */
105: PetscCall(PetscMalloc1(count, &fbuf));
106: PetscCall(PetscBinaryRead(fdes, fbuf, count, NULL, PETSC_FLOAT));
107: PetscCall(PetscByteSwap(fbuf, PETSC_FLOAT, count));
108: for (PetscInt i = 0; i < count; i++) ((PetscScalar *)data)[i] = fbuf[i];
109: PetscCall(PetscFree(fbuf));
110: } else {
111: PetscCall(PetscViewerASCIIRead(viewer, data, count, NULL, dtype));
112: }
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode DMPlexCreateFluent_ReadSection(PetscViewer viewer, FluentSection *s)
117: {
118: char buffer[PETSC_MAX_PATH_LEN];
119: int snum;
121: PetscFunctionBegin;
122: /* Fast-forward to next section and derive its index */
123: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
124: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ' '));
125: snum = sscanf(buffer, "%d", &s->index);
126: /* If we can't match an index return -1 to signal end-of-file */
127: if (snum < 1) {
128: s->index = -1;
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: if (s->index == 0) { /* Comment */
133: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
135: } else if (s->index == 2) { /* Dimension */
136: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
137: snum = sscanf(buffer, "%d", &s->nd);
138: PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
140: } else if (s->index == 10 || s->index == 2010) { /* Vertices */
141: PetscInt numClosingParens = 0;
143: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
144: snum = sscanf(buffer, "(%x %x %x %d %d)", &s->zoneID, &s->first, &s->last, &s->type, &s->nd);
145: PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
146: if (s->zoneID > 0) {
147: PetscInt numCoords = s->last - s->first + 1;
148: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
149: PetscCall(PetscMalloc1(s->nd * numCoords, (PetscScalar **)&s->data));
150: PetscCall(DMPlexCreateFluent_ReadValues(viewer, s->data, s->nd * numCoords, PETSC_SCALAR, s->index == 2010 ? PETSC_TRUE : PETSC_FALSE, &numClosingParens));
151: if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
152: else --numClosingParens;
153: }
154: if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
155: else --numClosingParens;
156: PetscCheck(!numClosingParens, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
157: } else if (s->index == 12 || s->index == 2012) { /* Cells */
158: PetscInt numClosingParens = 0;
160: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
161: snum = sscanf(buffer, "(%x", &s->zoneID);
162: PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
163: if (s->zoneID == 0) { /* Header section */
164: snum = sscanf(buffer, "(%x %x %x %d)", &s->zoneID, &s->first, &s->last, &s->nd);
165: PetscCheck(snum == 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
166: } else { /* Data section */
167: snum = sscanf(buffer, "(%x %x %x %d %d)", &s->zoneID, &s->first, &s->last, &s->type, &s->nd);
168: PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
169: if (s->nd == 0) {
170: /* Read cell type definitions for mixed cells */
171: PetscInt numCells = s->last - s->first + 1;
172: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
173: PetscCall(PetscMalloc1(numCells, (PetscInt **)&s->data));
174: PetscCall(DMPlexCreateFluent_ReadValues(viewer, s->data, numCells, PETSC_INT, s->index == 2012 ? PETSC_TRUE : PETSC_FALSE, &numClosingParens));
175: if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
176: else --numClosingParens;
177: }
178: }
179: if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
180: else --numClosingParens;
181: PetscCheck(!numClosingParens, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
182: } else if (s->index == 13 || s->index == 2013) { /* Faces */
183: PetscInt numClosingParens = 0;
185: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
186: snum = sscanf(buffer, "(%x", &s->zoneID);
187: PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
188: if (s->zoneID == 0) { /* Header section */
189: snum = sscanf(buffer, "(%x %x %x %d)", &s->zoneID, &s->first, &s->last, &s->nd);
190: PetscCheck(snum == 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
191: } else { /* Data section */
192: PetscInt numEntries, numFaces, maxsize = 0, offset = 0;
194: snum = sscanf(buffer, "(%x %x %x %d %d)", &s->zoneID, &s->first, &s->last, &s->type, &s->nd);
195: PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
196: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
197: switch (s->nd) {
198: case 0:
199: numEntries = PETSC_DETERMINE;
200: break;
201: case 2:
202: numEntries = 2 + 2;
203: break; /* linear */
204: case 3:
205: numEntries = 2 + 3;
206: break; /* triangular */
207: case 4:
208: numEntries = 2 + 4;
209: break; /* quadrilateral */
210: default:
211: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file");
212: }
213: numFaces = s->last - s->first + 1;
214: if (numEntries != PETSC_DETERMINE) {
215: /* Allocate space only if we already know the size of the block */
216: PetscCall(PetscMalloc1(numEntries * numFaces, (PetscInt **)&s->data));
217: }
218: for (PetscInt f = 0; f < numFaces; f++) {
219: if (s->nd == 0) {
220: /* Determine the size of the block for "mixed" facets */
221: PetscInt numFaceVert = 0;
222: PetscCall(DMPlexCreateFluent_ReadValues(viewer, &numFaceVert, 1, PETSC_INT, s->index == 2013 ? PETSC_TRUE : PETSC_FALSE, &numClosingParens));
223: if (!f) {
224: maxsize = (numFaceVert + 3) * numFaces;
225: PetscCall(PetscMalloc1(maxsize, (PetscInt **)&s->data));
226: } else {
227: if (offset + numFaceVert + 3 >= maxsize) {
228: PetscInt *tmp;
230: PetscCall(PetscMalloc1(maxsize * 2, &tmp));
231: PetscCall(PetscArraycpy(tmp, (PetscInt *)s->data, maxsize));
232: PetscCall(PetscFree(s->data));
233: maxsize *= 2;
234: s->data = tmp;
235: }
236: }
237: ((PetscInt *)s->data)[offset] = numFaceVert;
238: ++offset;
239: numEntries = numFaceVert + 2;
240: }
241: PetscCall(DMPlexCreateFluent_ReadValues(viewer, &(((PetscInt *)s->data)[offset]), numEntries, PETSC_INT, s->index == 2013 ? PETSC_TRUE : PETSC_FALSE, &numClosingParens));
242: offset += numEntries;
243: }
244: if (s->nd != 0) PetscCall(PetscMPIIntCast(numEntries - 2, &s->nd));
245: if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
246: else --numClosingParens;
247: }
248: if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
249: else --numClosingParens;
250: PetscCheck(!numClosingParens, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
251: } else if (s->index == 39) { /* Label information */
252: char labelName[PETSC_MAX_PATH_LEN];
253: char caseName[PETSC_MAX_PATH_LEN];
255: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
256: snum = sscanf(buffer, "(%u %s %s %d)", &s->zoneID, caseName, labelName, &s->nd);
257: PetscCheck(snum == 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file: %d", snum);
258: PetscInt depth = 1;
259: do {
260: /* Match parentheses when parsing unknown sections */
261: do PetscCall(PetscViewerRead(viewer, &buffer[0], 1, NULL, PETSC_CHAR));
262: while (buffer[0] != '(' && buffer[0] != ')');
263: if (buffer[0] == '(') depth++;
264: if (buffer[0] == ')') depth--;
265: } while (depth > 0);
266: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '\n'));
267: PetscCall(PetscStrallocpy(labelName, (char **)&s->data));
268: PetscCall(PetscInfo((PetscObject)viewer, "CASE: Zone ID %u is label %s\n", s->zoneID, labelName));
269: } else { /* Unknown section type */
270: PetscInt depth = 1;
271: do {
272: /* Match parentheses when parsing unknown sections */
273: do PetscCall(PetscViewerRead(viewer, &buffer[0], 1, NULL, PETSC_CHAR));
274: while (buffer[0] != '(' && buffer[0] != ')');
275: if (buffer[0] == '(') depth++;
276: if (buffer[0] == ')') depth--;
277: } while (depth > 0);
278: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '\n'));
279: }
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: // Inserts point `face` with orientation `ornt` into the cone of point `cell` at position `c`, which is the first empty slot
284: static PetscErrorCode InsertFace(DM dm, PetscInt cell, PetscInt face, PetscInt ornt)
285: {
286: const PetscInt *cone;
287: PetscInt coneSize, c;
289: PetscFunctionBegin;
290: PetscCall(DMPlexGetCone(dm, cell, &cone));
291: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
292: for (c = 0; c < coneSize; ++c)
293: if (cone[c] < 0) break;
294: PetscCheck(c < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be inserted in cone of cell %" PetscInt_FMT " with size %" PetscInt_FMT, face, cell, coneSize);
295: PetscCall(DMPlexInsertCone(dm, cell, c, face));
296: PetscCall(DMPlexInsertConeOrientation(dm, cell, c, ornt));
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: static PetscErrorCode ReorderPolygon(DM dm, PetscInt cell)
301: {
302: const PetscInt *cone, *ornt;
303: PetscInt coneSize, newCone[16], newOrnt[16];
305: PetscFunctionBegin;
306: PetscCall(DMPlexGetOrientedCone(dm, cell, &cone, &ornt));
307: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
308: newCone[0] = cone[0];
309: newOrnt[0] = ornt[0];
310: for (PetscInt c = 1; c < coneSize; ++c) {
311: const PetscInt *fcone;
312: PetscInt firstVertex, lastVertex, c2;
314: PetscCall(DMPlexGetCone(dm, newCone[c - 1], &fcone));
315: lastVertex = newOrnt[c - 1] ? fcone[0] : fcone[1];
316: for (c2 = 0; c2 < coneSize; ++c2) {
317: const PetscInt *fcone2;
319: PetscCall(DMPlexGetCone(dm, cone[c2], &fcone2));
320: firstVertex = ornt[c2] ? fcone2[1] : fcone2[0];
321: if (lastVertex == firstVertex) {
322: // Point `cell` matched point `lastVertex` on face `cone[c2]` with orientation `ornt[c2]`
323: break;
324: }
325: }
326: 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);
327: newCone[c] = cone[c2];
328: newOrnt[c] = ornt[c2];
329: }
330: {
331: const PetscInt *fcone, *fcone2;
332: PetscInt vertex, vertex2;
334: PetscCall(DMPlexGetCone(dm, newCone[coneSize - 1], &fcone));
335: PetscCall(DMPlexGetCone(dm, newCone[0], &fcone2));
336: vertex = newOrnt[coneSize - 1] ? fcone[0] : fcone[1];
337: vertex2 = newOrnt[0] ? fcone2[1] : fcone2[0];
338: PetscCheck(vertex == vertex2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " did not match at the endpoint", cell);
339: }
340: PetscCall(DMPlexSetCone(dm, cell, newCone));
341: PetscCall(DMPlexSetConeOrientation(dm, cell, newOrnt));
342: PetscCall(DMPlexRestoreOrientedCone(dm, cell, &cone, &ornt));
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: static PetscErrorCode ReorderTetrahedron(PetscViewer viewer, DM dm, PetscInt cell)
347: {
348: const PetscInt *cone, *ornt, *fcone, *fornt, *farr, faces[4] = {0, 1, 3, 2};
349: PetscInt newCone[16], newOrnt[16];
351: PetscFunctionBegin;
352: PetscCall(DMPlexGetOrientedCone(dm, cell, &cone, &ornt));
353: newCone[0] = cone[0];
354: newOrnt[0] = ornt[0];
355: PetscCall(DMPlexGetOrientedCone(dm, newCone[0], &fcone, &fornt));
356: farr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, newOrnt[0]);
357: // Loop over each edge in the initial triangle
358: for (PetscInt e = 0; e < 3; ++e) {
359: const PetscInt edge = fcone[farr[e * 2 + 0]], eornt = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[e * 2 + 1], fornt[farr[e * 2 + 0]]);
360: PetscInt c;
362: // Loop over each remaining face in the tetrahedron
363: // On face `newCone[0]`, trying to match edge `edge` with final orientation `eornt` to an edge on another face
364: for (c = 1; c < 4; ++c) {
365: const PetscInt *fcone2, *fornt2, *farr2;
366: PetscInt c2;
367: PetscBool flip = PETSC_FALSE;
369: // Checking face `cone[c]` with orientation `ornt[c]`
370: PetscCall(DMPlexGetOrientedCone(dm, cone[c], &fcone2, &fornt2));
371: farr2 = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, ornt[c]);
372: // Check for edge
373: for (c2 = 0; c2 < 3; ++c2) {
374: const PetscInt edge2 = fcone2[farr2[c2 * 2 + 0]], eornt2 = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr2[c2 * 2 + 1], fornt2[farr2[c2 * 2 + 0]]);
375: // Trying to match edge `edge2` with final orientation `eornt2`
376: if (edge == edge2) {
377: //PetscCheck(eornt == -(eornt2 + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell % " PetscInt_FMT " edge %" PetscInt_FMT " (%" PetscInt_FMT ") found twice with the same orientation in face %" PetscInt_FMT " edge %" PetscInt_FMT, cell, edge, e, c, c2);
378: // Matched face `newCone[0]` with orientation `newOrnt[0]` to face `cone[c]` with orientation `ornt[c]` along edge `edge`
379: PetscCall(PetscInfo((PetscObject)viewer, "CASE: Matched cell %" PetscInt_FMT " edge %" PetscInt_FMT "/%" PetscInt_FMT " (%" PetscInt_FMT ") to face %" PetscInt_FMT "/%" PetscInt_FMT " edge %" PetscInt_FMT " (%" PetscInt_FMT ")\n", cell, edge, e, eornt, cone[c], c, c2, eornt2));
380: flip = eornt != -(eornt2 + 1) ? PETSC_TRUE : PETSC_FALSE;
381: break;
382: }
383: }
384: if (c2 < 3) {
385: newCone[faces[e + 1]] = cone[c];
386: // Compute new orientation of face based on which edge was matched (only the first edge matches a side different from 0)
387: // Face 1 should match its edge 2
388: // Face 2 should match its edge 0
389: // Face 3 should match its edge 0
390: if (flip) {
391: newOrnt[faces[e + 1]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_TRIANGLE, -((c2 + (!e ? 1 : 2)) % 3 + 1), ornt[c]);
392: } else {
393: newOrnt[faces[e + 1]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_TRIANGLE, !e ? (c2 + 1) % 3 : c2, ornt[c]);
394: }
395: break;
396: }
397: }
398: 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);
399: }
400: PetscCall(DMPlexRestoreOrientedCone(dm, newCone[0], &fcone, &fornt));
401: PetscCall(DMPlexSetCone(dm, cell, newCone));
402: PetscCall(DMPlexSetConeOrientation(dm, cell, newOrnt));
403: PetscCall(DMPlexRestoreOrientedCone(dm, cell, &cone, &ornt));
404: PetscFunctionReturn(PETSC_SUCCESS);
405: }
407: static PetscErrorCode ReorderHexahedron(DM dm, PetscInt cell)
408: {
409: const PetscInt *cone, *ornt, *fcone, *fornt, *farr;
410: const PetscInt faces[6] = {0, 5, 3, 4, 2, 1};
411: PetscInt used[6] = {1, 0, 0, 0, 0, 0};
412: PetscInt newCone[16], newOrnt[16];
414: PetscFunctionBegin;
415: PetscCall(DMPlexGetOrientedCone(dm, cell, &cone, &ornt));
416: newCone[0] = cone[0];
417: newOrnt[0] = ornt[0];
418: PetscCall(DMPlexGetOrientedCone(dm, newCone[0], &fcone, &fornt));
419: farr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, newOrnt[0]);
420: // Loop over each edge in the initial quadrilateral
421: for (PetscInt e = 0; e < 4; ++e) {
422: const PetscInt edge = fcone[farr[e * 2 + 0]], eornt = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[e * 2 + 1], fornt[farr[e * 2 + 0]]);
423: PetscInt c;
425: // Loop over each remaining face in the hexahedron
426: // On face `newCone[0]`, trying to match edge `edge` with final orientation `eornt` to an edge on another face
427: for (c = 1; c < 6; ++c) {
428: const PetscInt *fcone2, *fornt2, *farr2;
429: PetscInt c2;
431: // Checking face `cone[c]` with orientation `ornt[c]`
432: PetscCall(DMPlexGetOrientedCone(dm, cone[c], &fcone2, &fornt2));
433: farr2 = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, ornt[c]);
434: // Check for edge
435: for (c2 = 0; c2 < 4; ++c2) {
436: const PetscInt edge2 = fcone2[farr2[c2 * 2 + 0]], eornt2 = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr2[c2 * 2 + 1], fornt2[farr2[c2 * 2 + 0]]);
437: // Trying to match edge `edge2` with final orientation `eornt2`
438: if (edge == edge2) {
439: PetscCheck(eornt == -(eornt2 + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " found twice with the same orientation", edge);
440: // Matched face `newCone[0]` with orientation `newOrnt[0]` to face `cone[c]` with orientation `ornt[c]` along edge `edge`
441: break;
442: }
443: }
444: if (c2 < 4) {
445: used[c] = 1;
446: newCone[faces[e + 1]] = cone[c];
447: // Compute new orientation of face based on which edge was matched (only the first edge matches a side different from 0)
448: newOrnt[faces[e + 1]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_QUADRILATERAL, !e ? (c2 + 1) % 4 : c2, ornt[c]);
449: break;
450: }
451: }
452: 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);
453: }
454: PetscCall(DMPlexRestoreOrientedCone(dm, newCone[0], &fcone, &fornt));
455: // Add last face
456: {
457: PetscInt c, c2;
459: for (c = 1; c < 6; ++c)
460: if (!used[c]) break;
461: PetscCheck(c < 6, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " could not find an available face", cell);
462: // Match first edge to 3rd edge in newCone[2]
463: {
464: const PetscInt *fcone2, *fornt2, *farr2;
466: PetscCall(DMPlexGetOrientedCone(dm, newCone[2], &fcone, &fornt));
467: farr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, newOrnt[2]);
468: PetscCall(DMPlexGetOrientedCone(dm, cone[c], &fcone2, &fornt2));
469: farr2 = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, ornt[c]);
471: const PetscInt e = 2;
472: const PetscInt edge = fcone[farr[e * 2 + 0]], eornt = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[e * 2 + 1], fornt[farr[e * 2 + 0]]);
473: // Trying to match edge `edge` with final orientation `eornt` of face `newCone[2]` to some edge of face `cone[c]` with orientation `ornt[c]`
474: for (c2 = 0; c2 < 4; ++c2) {
475: const PetscInt edge2 = fcone2[farr2[c2 * 2 + 0]], eornt2 = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr2[c2 * 2 + 1], fornt2[farr2[c2 * 2 + 0]]);
476: // Trying to match edge `edge2` with final orientation `eornt2`
477: if (edge == edge2) {
478: PetscCheck(eornt == -(eornt2 + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " found twice with the same orientation", edge);
479: // Matched face `newCone[2]` with orientation `newOrnt[2]` to face `cone[c]` with orientation `ornt[c]` along edge `edge`
480: break;
481: }
482: }
483: PetscCheck(c2 < 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not fit last face in");
484: }
485: newCone[faces[5]] = cone[c];
486: // Compute new orientation of face based on which edge was matched
487: newOrnt[faces[5]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_QUADRILATERAL, c2, ornt[c]);
488: PetscCall(DMPlexRestoreOrientedCone(dm, newCone[0], &fcone, &fornt));
489: }
490: PetscCall(DMPlexSetCone(dm, cell, newCone));
491: PetscCall(DMPlexSetConeOrientation(dm, cell, newOrnt));
492: PetscCall(DMPlexRestoreOrientedCone(dm, cell, &cone, &ornt));
493: PetscFunctionReturn(PETSC_SUCCESS);
494: }
496: // {0, 1, 2}, {3, 4, 5}, {0, 2, 4, 3}, {2, 1, 5, 4}, {1, 0, 3, 5}
497: static PetscErrorCode ReorderWedge(DM dm, PetscInt cell)
498: {
499: const PetscInt *cone, *ornt, *fcone, *fornt, *farr;
500: const PetscInt faces[5] = {0, 4, 3, 2, 1};
501: PetscInt used[5] = {0, 0, 0, 0, 0};
502: PetscInt newCone[16], newOrnt[16], cS, bottom = 0;
504: PetscFunctionBegin;
505: PetscCall(DMPlexGetConeSize(dm, cell, &cS));
506: PetscCall(DMPlexGetOrientedCone(dm, cell, &cone, &ornt));
507: for (PetscInt c = 0; c < cS; ++c) {
508: DMPolytopeType ct;
510: PetscCall(DMPlexGetCellType(dm, cone[c], &ct));
511: if (ct == DM_POLYTOPE_TRIANGLE) {
512: bottom = c;
513: break;
514: }
515: }
516: used[bottom] = 1;
517: newCone[0] = cone[bottom];
518: newOrnt[0] = ornt[bottom];
519: PetscCall(DMPlexGetOrientedCone(dm, newCone[0], &fcone, &fornt));
520: farr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, newOrnt[0]);
521: // Loop over each edge in the initial triangle
522: for (PetscInt e = 0; e < 3; ++e) {
523: const PetscInt edge = fcone[farr[e * 2 + 0]], eornt = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[e * 2 + 1], fornt[farr[e * 2 + 0]]);
524: PetscInt c;
526: // Loop over each remaining face in the prism
527: // On face `newCone[0]`, trying to match edge `edge` with final orientation `eornt` to an edge on another face
528: for (c = 0; c < 5; ++c) {
529: const PetscInt *fcone2, *fornt2, *farr2;
530: DMPolytopeType ct;
531: PetscInt c2;
533: if (c == bottom) continue;
534: PetscCall(DMPlexGetCellType(dm, cone[c], &ct));
535: if (ct != DM_POLYTOPE_QUADRILATERAL) continue;
536: // Checking face `cone[c]` with orientation `ornt[c]`
537: PetscCall(DMPlexGetOrientedCone(dm, cone[c], &fcone2, &fornt2));
538: farr2 = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, ornt[c]);
539: // Check for edge
540: for (c2 = 0; c2 < 4; ++c2) {
541: const PetscInt edge2 = fcone2[farr2[c2 * 2 + 0]], eornt2 = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr2[c2 * 2 + 1], fornt2[farr2[c2 * 2 + 0]]);
542: // Trying to match edge `edge2` with final orientation `eornt2`
543: if (edge == edge2) {
544: PetscCheck(eornt == -(eornt2 + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " found twice with the same orientation", edge);
545: // Matched face `newCone[0]` with orientation `newOrnt[0]` to face `cone[c]` with orientation `ornt[c]` along edge `edge`
546: break;
547: }
548: }
549: if (c2 < 4) {
550: used[c] = 1;
551: newCone[faces[e + 1]] = cone[c];
552: // Compute new orientation of face based on which edge was matched, edge 0 should always match the bottom
553: newOrnt[faces[e + 1]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_QUADRILATERAL, c2, ornt[c]);
554: break;
555: }
556: }
557: PetscCheck(c < 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " could not find a face match for edge %" PetscInt_FMT, cell, e);
558: }
559: PetscCall(DMPlexRestoreOrientedCone(dm, newCone[0], &fcone, &fornt));
560: // Add last face
561: {
562: PetscInt c, c2;
564: for (c = 0; c < 5; ++c)
565: if (!used[c]) break;
566: PetscCheck(c < 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " could not find an available face", cell);
567: // Match first edge to 3rd edge in newCone[2]
568: {
569: const PetscInt *fcone2, *fornt2, *farr2;
571: PetscCall(DMPlexGetOrientedCone(dm, newCone[2], &fcone, &fornt));
572: farr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, newOrnt[2]);
573: PetscCall(DMPlexGetOrientedCone(dm, cone[c], &fcone2, &fornt2));
574: farr2 = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, ornt[c]);
576: const PetscInt e = 2;
577: const PetscInt edge = fcone[farr[e * 2 + 0]], eornt = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[e * 2 + 1], fornt[farr[e * 2 + 0]]);
578: // Trying to match edge `edge` with final orientation `eornt` of face `newCone[2]` to some edge of face `cone[c]` with orientation `ornt[c]`
579: for (c2 = 0; c2 < 3; ++c2) {
580: const PetscInt edge2 = fcone2[farr2[c2 * 2 + 0]], eornt2 = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr2[c2 * 2 + 1], fornt2[farr2[c2 * 2 + 0]]);
581: // Trying to match edge `edge2` with final orientation `eornt2`
582: if (edge == edge2) {
583: PetscCheck(eornt == -(eornt2 + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " found twice with the same orientation", edge);
584: // Matched face `newCone[2]` with orientation `newOrnt[2]` to face `cone[c]` with orientation `ornt[c]` along edge `edge`
585: break;
586: }
587: }
588: PetscCheck(c2 < 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not fit last face in");
589: }
590: newCone[faces[4]] = cone[c];
591: // Compute new orientation of face based on which edge was matched
592: newOrnt[faces[4]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_TRIANGLE, c2, ornt[c]);
593: PetscCall(DMPlexRestoreOrientedCone(dm, newCone[0], &fcone, &fornt));
594: }
595: PetscCall(DMPlexSetCone(dm, cell, newCone));
596: PetscCall(DMPlexSetConeOrientation(dm, cell, newOrnt));
597: PetscCall(DMPlexRestoreOrientedCone(dm, cell, &cone, &ornt));
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: static PetscErrorCode ReorderCell(PetscViewer viewer, DM dm, PetscInt cell, DMPolytopeType ct)
602: {
603: PetscFunctionBegin;
604: switch (ct) {
605: case DM_POLYTOPE_TRIANGLE:
606: case DM_POLYTOPE_QUADRILATERAL:
607: PetscCall(ReorderPolygon(dm, cell));
608: break;
609: case DM_POLYTOPE_TETRAHEDRON:
610: PetscCall(ReorderTetrahedron(viewer, dm, cell));
611: break;
612: case DM_POLYTOPE_HEXAHEDRON:
613: PetscCall(ReorderHexahedron(dm, cell));
614: break;
615: case DM_POLYTOPE_TRI_PRISM:
616: PetscCall(ReorderWedge(dm, cell));
617: break;
618: default:
619: PetscCheck(0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Celltype %s is unsupported", DMPolytopeTypes[ct]);
620: break;
621: }
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
625: static PetscErrorCode GetNumCellFaces(int nd, PetscInt *numCellFaces, DMPolytopeType *ct)
626: {
627: PetscFunctionBegin;
628: *ct = DM_POLYTOPE_POINT;
629: switch (nd) {
630: case 0:
631: *numCellFaces = PETSC_DETERMINE;
632: break;
633: case 1:
634: *numCellFaces = 3;
635: *ct = DM_POLYTOPE_TRIANGLE;
636: break;
637: case 2:
638: *numCellFaces = 4;
639: *ct = DM_POLYTOPE_TETRAHEDRON;
640: break;
641: case 3:
642: *numCellFaces = 4;
643: *ct = DM_POLYTOPE_QUADRILATERAL;
644: break;
645: case 4:
646: *numCellFaces = 6;
647: *ct = DM_POLYTOPE_HEXAHEDRON;
648: break;
649: case 5:
650: *numCellFaces = 5;
651: *ct = DM_POLYTOPE_PYRAMID;
652: break;
653: case 6:
654: *numCellFaces = 5;
655: *ct = DM_POLYTOPE_TRI_PRISM;
656: break;
657: default:
658: *numCellFaces = PETSC_DETERMINE;
659: }
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: /*@C
664: DMPlexCreateFluent - Create a `DMPLEX` mesh from a Fluent mesh file <http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm>.
666: Collective
668: Input Parameters:
669: + comm - The MPI communicator
670: . viewer - The `PetscViewer` associated with a Fluent mesh file
671: - interpolate - Create faces and edges in the mesh
673: Output Parameter:
674: . dm - The `DM` object representing the mesh
676: Level: beginner
678: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`
679: @*/
680: PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
681: {
682: PetscInt dim = PETSC_DETERMINE;
683: PetscInt numCells = 0;
684: PetscInt numVertices = 0;
685: PetscInt *cellSizes = NULL;
686: DMPolytopeType *cellTypes = NULL;
687: PetscInt numFaces = 0;
688: PetscInt *faces = NULL;
689: PetscInt *faceSizes = NULL;
690: PetscInt *faceAdjCell = NULL;
691: PetscInt *cellVertices = NULL;
692: unsigned int *faceZoneIDs = NULL;
693: DMLabel faceSets = NULL;
694: DMLabel *zoneLabels = NULL;
695: const char **zoneNames = NULL;
696: unsigned int maxZoneID = 0;
697: PetscScalar *coordsIn = NULL;
698: PetscScalar *coords;
699: PetscSection coordSection;
700: Vec coordinates;
701: PetscInt coordSize, maxFaceSize = 0, totFaceVert = 0, f;
702: PetscMPIInt rank;
704: PetscFunctionBegin;
705: PetscCallMPI(MPI_Comm_rank(comm, &rank));
707: if (rank == 0) {
708: FluentSection s;
710: s.data = NULL;
711: numFaces = PETSC_DETERMINE;
712: do {
713: PetscCall(DMPlexCreateFluent_ReadSection(viewer, &s));
714: if (s.index == 2) { /* Dimension */
715: dim = s.nd;
716: PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found dimension: %" PetscInt_FMT "\n", dim));
717: } else if (s.index == 10 || s.index == 2010) { /* Vertices */
718: if (s.zoneID == 0) {
719: numVertices = s.last;
720: PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found number of vertices: %" PetscInt_FMT "\n", numVertices));
721: } else {
722: PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found vertex coordinates\n"));
723: PetscCheck(!coordsIn, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files");
724: coordsIn = (PetscScalar *)s.data;
725: }
727: } else if (s.index == 12 || s.index == 2012) { /* Cells */
728: if (s.zoneID == 0) {
729: numCells = s.last;
730: PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found number of cells %" PetscInt_FMT "\n", numCells));
731: } else {
732: PetscCall(PetscMalloc2(numCells, &cellSizes, numCells, &cellTypes));
733: for (PetscInt c = 0; c < numCells; ++c) PetscCall(GetNumCellFaces(s.nd ? s.nd : (int)((PetscInt *)s.data)[c], &cellSizes[c], &cellTypes[c]));
734: PetscCall(PetscFree(s.data));
735: PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found number of cell faces %" PetscInt_FMT "\n", numCells && s.nd ? cellSizes[0] : 0));
736: }
737: } else if (s.index == 13 || s.index == 2013) { /* Facets */
738: if (s.zoneID == 0) { /* Header section */
739: numFaces = (PetscInt)(s.last - s.first + 1);
740: PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found number of faces %" PetscInt_FMT " face vertices: %d\n", numFaces, s.nd));
741: } else { /* Data section */
742: PetscInt *tmp;
743: PetscInt totSize = 0, offset = 0, doffset;
745: PetscCheck(numFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file");
746: if (!faceZoneIDs) PetscCall(PetscMalloc3(numFaces, &faceSizes, numFaces * 2, &faceAdjCell, numFaces, &faceZoneIDs));
747: // Record the zoneID and face size for each face set
748: for (unsigned int z = s.first - 1; z < s.last; z++) {
749: faceZoneIDs[z] = s.zoneID;
750: if (s.nd) {
751: faceSizes[z] = s.nd;
752: } else {
753: faceSizes[z] = ((PetscInt *)s.data)[offset];
754: offset += faceSizes[z] + 3;
755: }
756: totSize += faceSizes[z];
757: maxFaceSize = PetscMax(maxFaceSize, faceSizes[z]);
758: }
760: offset = totFaceVert;
761: doffset = s.nd ? 0 : 1;
762: PetscCall(PetscMalloc1(totFaceVert + totSize, &tmp));
763: if (faces) PetscCall(PetscArraycpy(tmp, faces, totFaceVert));
764: PetscCall(PetscFree(faces));
765: totFaceVert += totSize;
766: faces = tmp;
768: // Record face vertices and adjacent faces
769: const PetscInt Nfz = s.last - s.first + 1;
770: for (PetscInt f = 0; f < Nfz; ++f) {
771: const PetscInt face = f + s.first - 1;
772: const PetscInt faceSize = faceSizes[face];
774: for (PetscInt v = 0; v < faceSize; ++v) faces[offset + v] = ((PetscInt *)s.data)[doffset + v];
775: faceAdjCell[face * 2 + 0] = ((PetscInt *)s.data)[doffset + faceSize + 0];
776: faceAdjCell[face * 2 + 1] = ((PetscInt *)s.data)[doffset + faceSize + 1];
777: offset += faceSize;
778: doffset += faceSize + (s.nd ? 2 : 3);
779: }
780: PetscCall(PetscFree(s.data));
781: }
782: } else if (s.index == 39) { /* Label information */
783: if (s.zoneID >= maxZoneID) {
784: DMLabel *tmpL;
785: const char **tmp;
786: unsigned int newmax = maxZoneID + 1;
788: while (newmax < s.zoneID + 1) newmax *= 2;
789: PetscCall(PetscCalloc2(newmax, &tmp, newmax, &tmpL));
790: for (PetscInt i = 0; i < (PetscInt)maxZoneID; ++i) {
791: tmp[i] = zoneNames[i];
792: tmpL[i] = zoneLabels[i];
793: }
794: maxZoneID = newmax;
795: PetscCall(PetscFree2(zoneNames, zoneLabels));
796: zoneNames = tmp;
797: zoneLabels = tmpL;
798: }
799: zoneNames[s.zoneID] = (const char *)s.data;
800: }
801: } while (s.index >= 0);
802: }
803: PetscCallMPI(MPI_Bcast(&dim, 1, MPIU_INT, 0, comm));
804: PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension");
806: /* Allocate cell-vertex mesh */
807: PetscCall(DMCreate(comm, dm));
808: PetscCall(DMSetType(*dm, DMPLEX));
809: PetscCall(DMSetDimension(*dm, dim));
810: // We do not want this label automatically computed, instead we fill it here
811: PetscCall(DMCreateLabel(*dm, "celltype"));
812: PetscCall(DMPlexSetChart(*dm, 0, numCells + numFaces + numVertices));
813: if (rank == 0) {
814: PetscCheck(numCells >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown number of cells in Fluent file");
815: for (PetscInt c = 0; c < numCells; ++c) {
816: PetscCall(DMPlexSetConeSize(*dm, c, cellSizes[c]));
817: PetscCall(DMPlexSetCellType(*dm, c, cellTypes[c]));
818: }
819: for (PetscInt v = numCells; v < numCells + numVertices; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT));
820: for (PetscInt f = 0; f < numFaces; ++f) {
821: DMPolytopeType ct;
823: switch (faceSizes[f]) {
824: case 2:
825: ct = DM_POLYTOPE_SEGMENT;
826: break;
827: case 3:
828: ct = DM_POLYTOPE_TRIANGLE;
829: break;
830: case 4:
831: ct = DM_POLYTOPE_QUADRILATERAL;
832: break;
833: default:
834: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file with cone size %" PetscInt_FMT, faceSizes[f]);
835: }
836: PetscCall(DMPlexSetConeSize(*dm, f + numCells + numVertices, faceSizes[f]));
837: PetscCall(DMPlexSetCellType(*dm, f + numCells + numVertices, ct));
838: }
839: }
840: PetscCall(DMSetUp(*dm));
842: if (rank == 0 && faces) {
843: PetscSection s;
844: PetscInt *cones, csize, foffset = 0;
846: PetscCall(DMPlexGetCones(*dm, &cones));
847: PetscCall(DMPlexGetConeSection(*dm, &s));
848: PetscCall(PetscSectionGetConstrainedStorageSize(s, &csize));
849: for (PetscInt c = 0; c < csize; ++c) cones[c] = -1;
850: for (PetscInt f = 0; f < numFaces; f++) {
851: const PetscInt cl = faceAdjCell[f * 2 + 0] - 1;
852: const PetscInt cr = faceAdjCell[f * 2 + 1] - 1;
853: const PetscInt face = f + numCells + numVertices;
854: PetscInt fcone[16];
856: // How could Fluent define the outward normal differently? Is there no end to the pain?
857: if (dim == 3) {
858: if (cl >= 0) PetscCall(InsertFace(*dm, cl, face, -1));
859: if (cr >= 0) PetscCall(InsertFace(*dm, cr, face, 0));
860: } else {
861: if (cl >= 0) PetscCall(InsertFace(*dm, cl, face, 0));
862: if (cr >= 0) PetscCall(InsertFace(*dm, cr, face, -1));
863: }
864: PetscCheck(faceSizes[f] < 16, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of face vertices %" PetscInt_FMT " exceeds temporary storage", faceSizes[f]);
865: for (PetscInt v = 0; v < faceSizes[f]; ++v) fcone[v] = faces[foffset + v] + numCells - 1;
866: foffset += faceSizes[f];
867: PetscCall(DMPlexSetCone(*dm, face, fcone));
868: }
869: }
870: PetscCall(DMPlexSymmetrize(*dm));
871: PetscCall(DMPlexStratify(*dm));
872: if (dim == 3) {
873: DM idm;
875: PetscCall(DMCreate(PetscObjectComm((PetscObject)*dm), &idm));
876: PetscCall(DMSetType(idm, DMPLEX));
877: PetscCall(DMSetDimension(idm, dim));
878: PetscCall(DMPlexInterpolateFaces_Internal(*dm, 1, idm));
879: PetscCall(DMDestroy(dm));
880: *dm = idm;
881: }
882: PetscCall(DMViewFromOptions(*dm, NULL, "-cas_dm_view"));
883: if (rank == 0 && faces) {
884: for (PetscInt c = 0; c < numCells; ++c) PetscCall(ReorderCell(viewer, *dm, c, cellTypes[c]));
885: }
887: if (rank == 0 && faces) {
888: PetscInt joinSize, meetSize, *fverts, cells[2];
889: const PetscInt *join, *meet;
890: PetscInt foffset = 0;
892: PetscCall(PetscMalloc1(maxFaceSize, &fverts));
893: /* Mark facets by finding the full join of all adjacent vertices */
894: for (f = 0; f < numFaces; f++) {
895: const PetscInt cl = faceAdjCell[f * 2 + 0] - 1;
896: const PetscInt cr = faceAdjCell[f * 2 + 1] - 1;
897: const PetscInt id = (PetscInt)faceZoneIDs[f];
899: if (cl > 0 && cr > 0) {
900: /* If we know both adjoining cells we can use a single-level meet */
901: cells[0] = cl;
902: cells[1] = cr;
903: PetscCall(DMPlexGetMeet(*dm, 2, cells, &meetSize, &meet));
904: 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);
905: PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", meet[0], id));
906: if (zoneNames && zoneNames[id]) PetscCall(DMSetLabelValue_Fast(*dm, &zoneLabels[id], zoneNames[id], meet[0], 1));
907: PetscCall(DMPlexRestoreMeet(*dm, meetSize, fverts, &meetSize, &meet));
908: } else {
909: for (PetscInt fi = 0; fi < faceSizes[f]; fi++) fverts[fi] = faces[foffset + fi] + numCells - 1;
910: PetscCall(DMPlexGetFullJoin(*dm, faceSizes[f], fverts, &joinSize, &join));
911: PetscCheck(joinSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %" PetscInt_FMT, f);
912: PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], id));
913: if (zoneNames && zoneNames[id]) PetscCall(DMSetLabelValue_Fast(*dm, &zoneLabels[id], zoneNames[id], join[0], 1));
914: PetscCall(DMPlexRestoreJoin(*dm, joinSize, fverts, &joinSize, &join));
915: }
916: foffset += faceSizes[f];
917: }
918: PetscCall(PetscFree(fverts));
919: }
921: { /* Create Face Sets label at all processes */
922: enum {
923: n = 1
924: };
925: PetscBool flag[n];
927: flag[0] = faceSets ? PETSC_TRUE : PETSC_FALSE;
928: PetscCallMPI(MPI_Bcast(flag, n, MPI_C_BOOL, 0, comm));
929: if (flag[0]) PetscCall(DMCreateLabel(*dm, "Face Sets"));
930: // TODO Code to create all the zone labels on each process
931: }
933: if (!interpolate) {
934: DM udm;
936: PetscCall(DMPlexUninterpolate(*dm, &udm));
937: PetscCall(DMDestroy(dm));
938: *dm = udm;
939: }
941: /* Read coordinates */
942: PetscCall(DMGetCoordinateSection(*dm, &coordSection));
943: PetscCall(PetscSectionSetNumFields(coordSection, 1));
944: PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dim));
945: PetscCall(PetscSectionSetChart(coordSection, numCells, numCells + numVertices));
946: for (PetscInt v = numCells; v < numCells + numVertices; ++v) {
947: PetscCall(PetscSectionSetDof(coordSection, v, dim));
948: PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dim));
949: }
950: PetscCall(PetscSectionSetUp(coordSection));
951: PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize));
952: PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates));
953: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
954: PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE));
955: PetscCall(VecSetType(coordinates, VECSTANDARD));
956: PetscCall(VecGetArray(coordinates, &coords));
957: if (rank == 0 && coordsIn) {
958: for (PetscInt v = 0; v < numVertices; ++v) {
959: for (PetscInt d = 0; d < dim; ++d) coords[v * dim + d] = coordsIn[v * dim + d];
960: }
961: }
962: PetscCall(VecRestoreArray(coordinates, &coords));
963: PetscCall(DMSetCoordinatesLocal(*dm, coordinates));
964: PetscCall(VecDestroy(&coordinates));
966: if (rank == 0) {
967: PetscCall(PetscFree(cellVertices));
968: PetscCall(PetscFree2(cellSizes, cellTypes));
969: PetscCall(PetscFree(faces));
970: PetscCall(PetscFree3(faceSizes, faceAdjCell, faceZoneIDs));
971: PetscCall(PetscFree(coordsIn));
972: if (zoneNames)
973: for (PetscInt i = 0; i < (PetscInt)maxZoneID; ++i) PetscCall(PetscFree(zoneNames[i]));
974: PetscCall(PetscFree2(zoneNames, zoneLabels));
975: }
976: PetscFunctionReturn(PETSC_SUCCESS);
977: }