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