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