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: /*@C
6: DMPlexCreateFluentFromFile - Create a `DMPLEX` mesh from a Fluent mesh file
8: Collective
10: Input Parameters:
11: + comm - The MPI communicator
12: . filename - Name of the Fluent mesh file
13: - interpolate - Create faces and edges in the mesh
15: Output Parameter:
16: . dm - The `DM` object representing the mesh
18: Level: beginner
20: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateFromFile()`, `DMPlexCreateFluent()`, `DMPlexCreate()`
21: @*/
22: PetscErrorCode DMPlexCreateFluentFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
23: {
24: PetscViewer viewer;
26: PetscFunctionBegin;
27: /* Create file viewer and build plex */
28: PetscCall(PetscViewerCreate(comm, &viewer));
29: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII));
30: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ));
31: PetscCall(PetscViewerFileSetName(viewer, filename));
32: PetscCall(DMPlexCreateFluent(comm, viewer, interpolate, dm));
33: PetscCall(PetscViewerDestroy(&viewer));
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: static PetscErrorCode DMPlexCreateFluent_ReadString(PetscViewer viewer, char *buffer, char delim)
38: {
39: PetscInt ret, i = 0;
41: PetscFunctionBegin;
42: do PetscCall(PetscViewerRead(viewer, &(buffer[i++]), 1, &ret, PETSC_CHAR));
43: while (ret > 0 && buffer[i - 1] != '\0' && buffer[i - 1] != delim);
44: if (!ret) buffer[i - 1] = '\0';
45: else buffer[i] = '\0';
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: static PetscErrorCode DMPlexCreateFluent_ReadValues(PetscViewer viewer, void *data, PetscInt count, PetscDataType dtype, PetscBool binary)
50: {
51: int fdes = 0;
52: FILE *file;
53: PetscInt i;
55: PetscFunctionBegin;
56: if (binary) {
57: /* Extract raw file descriptor to read binary block */
58: PetscCall(PetscViewerASCIIGetPointer(viewer, &file));
59: PetscCall(PetscFFlush(file));
60: fdes = fileno(file);
61: }
63: if (!binary && dtype == PETSC_INT) {
64: char cbuf[256];
65: unsigned int ibuf;
66: int snum;
67: /* Parse hexadecimal ascii integers */
68: for (i = 0; i < count; i++) {
69: PetscCall(PetscViewerRead(viewer, cbuf, 1, NULL, PETSC_STRING));
70: snum = sscanf(cbuf, "%x", &ibuf);
71: PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
72: ((PetscInt *)data)[i] = (PetscInt)ibuf;
73: }
74: } else if (binary && dtype == PETSC_INT) {
75: /* Always read 32-bit ints and cast to PetscInt */
76: int *ibuf;
77: PetscCall(PetscMalloc1(count, &ibuf));
78: PetscCall(PetscBinaryRead(fdes, ibuf, count, NULL, PETSC_ENUM));
79: PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, count));
80: for (i = 0; i < count; i++) ((PetscInt *)data)[i] = (PetscInt)(ibuf[i]);
81: PetscCall(PetscFree(ibuf));
83: } else if (binary && dtype == PETSC_SCALAR) {
84: float *fbuf;
85: /* Always read 32-bit floats and cast to PetscScalar */
86: PetscCall(PetscMalloc1(count, &fbuf));
87: PetscCall(PetscBinaryRead(fdes, fbuf, count, NULL, PETSC_FLOAT));
88: PetscCall(PetscByteSwap(fbuf, PETSC_FLOAT, count));
89: for (i = 0; i < count; i++) ((PetscScalar *)data)[i] = (PetscScalar)(fbuf[i]);
90: PetscCall(PetscFree(fbuf));
91: } else {
92: PetscCall(PetscViewerASCIIRead(viewer, data, count, NULL, dtype));
93: }
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: static PetscErrorCode DMPlexCreateFluent_ReadSection(PetscViewer viewer, FluentSection *s)
98: {
99: char buffer[PETSC_MAX_PATH_LEN];
100: int snum;
102: PetscFunctionBegin;
103: /* Fast-forward to next section and derive its index */
104: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
105: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ' '));
106: snum = sscanf(buffer, "%d", &(s->index));
107: /* If we can't match an index return -1 to signal end-of-file */
108: if (snum < 1) {
109: s->index = -1;
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: if (s->index == 0) { /* Comment */
114: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
116: } else if (s->index == 2) { /* Dimension */
117: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
118: snum = sscanf(buffer, "%d", &(s->nd));
119: PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
121: } else if (s->index == 10 || s->index == 2010) { /* Vertices */
122: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
123: snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
124: PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
125: if (s->zoneID > 0) {
126: PetscInt numCoords = s->last - s->first + 1;
127: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
128: PetscCall(PetscMalloc1(s->nd * numCoords, (PetscScalar **)&s->data));
129: PetscCall(DMPlexCreateFluent_ReadValues(viewer, s->data, s->nd * numCoords, PETSC_SCALAR, s->index == 2010 ? PETSC_TRUE : PETSC_FALSE));
130: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
131: }
132: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
134: } else if (s->index == 12 || s->index == 2012) { /* Cells */
135: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
136: snum = sscanf(buffer, "(%x", &(s->zoneID));
137: PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
138: if (s->zoneID == 0) { /* Header section */
139: snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
140: PetscCheck(snum == 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
141: } else { /* Data section */
142: snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
143: PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
144: if (s->nd == 0) {
145: /* Read cell type definitions for mixed cells */
146: PetscInt numCells = s->last - s->first + 1;
147: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
148: PetscCall(PetscMalloc1(numCells, (PetscInt **)&s->data));
149: PetscCall(DMPlexCreateFluent_ReadValues(viewer, s->data, numCells, PETSC_INT, s->index == 2012 ? PETSC_TRUE : PETSC_FALSE));
150: PetscCall(PetscFree(s->data));
151: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
152: }
153: }
154: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
156: } else if (s->index == 13 || s->index == 2013) { /* Faces */
157: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
158: snum = sscanf(buffer, "(%x", &(s->zoneID));
159: PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
160: if (s->zoneID == 0) { /* Header section */
161: snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
162: PetscCheck(snum == 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
163: } else { /* Data section */
164: PetscInt f, numEntries, numFaces;
165: snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
166: PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
167: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
168: switch (s->nd) {
169: case 0:
170: numEntries = PETSC_DETERMINE;
171: break;
172: case 2:
173: numEntries = 2 + 2;
174: break; /* linear */
175: case 3:
176: numEntries = 2 + 3;
177: break; /* triangular */
178: case 4:
179: numEntries = 2 + 4;
180: break; /* quadrilateral */
181: default:
182: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file");
183: }
184: numFaces = s->last - s->first + 1;
185: if (numEntries != PETSC_DETERMINE) {
186: /* Allocate space only if we already know the size of the block */
187: PetscCall(PetscMalloc1(numEntries * numFaces, (PetscInt **)&s->data));
188: }
189: for (f = 0; f < numFaces; f++) {
190: if (s->nd == 0) {
191: /* Determine the size of the block for "mixed" facets */
192: PetscInt numFaceVert = 0;
193: PetscCall(DMPlexCreateFluent_ReadValues(viewer, &numFaceVert, 1, PETSC_INT, s->index == 2013 ? PETSC_TRUE : PETSC_FALSE));
194: if (numEntries == PETSC_DETERMINE) {
195: numEntries = numFaceVert + 2;
196: PetscCall(PetscMalloc1(numEntries * numFaces, (PetscInt **)&s->data));
197: } else {
198: PetscCheck(numEntries == numFaceVert + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for mixed faces in Fluent files");
199: }
200: }
201: PetscCall(DMPlexCreateFluent_ReadValues(viewer, &(((PetscInt *)s->data)[f * numEntries]), numEntries, PETSC_INT, s->index == 2013 ? PETSC_TRUE : PETSC_FALSE));
202: }
203: s->nd = numEntries - 2;
204: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
205: }
206: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
208: } else { /* Unknown section type */
209: PetscInt depth = 1;
210: do {
211: /* Match parentheses when parsing unknown sections */
212: do PetscCall(PetscViewerRead(viewer, &(buffer[0]), 1, NULL, PETSC_CHAR));
213: while (buffer[0] != '(' && buffer[0] != ')');
214: if (buffer[0] == '(') depth++;
215: if (buffer[0] == ')') depth--;
216: } while (depth > 0);
217: PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '\n'));
218: }
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: /*@C
223: DMPlexCreateFluent - Create a `DMPLEX` mesh from a Fluent mesh file.
225: Collective
227: Input Parameters:
228: + comm - The MPI communicator
229: . viewer - The `PetscViewer` associated with a Fluent mesh file
230: - interpolate - Create faces and edges in the mesh
232: Output Parameter:
233: . dm - The `DM` object representing the mesh
235: Note:
236: http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm
238: Level: beginner
240: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`
241: @*/
242: PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
243: {
244: PetscMPIInt rank;
245: PetscInt c, v, dim = PETSC_DETERMINE, numCells = 0, numVertices = 0, numCellVertices = PETSC_DETERMINE;
246: PetscInt numFaces = PETSC_DETERMINE, f, numFaceEntries = PETSC_DETERMINE, numFaceVertices = PETSC_DETERMINE;
247: PetscInt *faces = NULL, *cellVertices = NULL, *faceZoneIDs = NULL;
248: DMLabel faceSets = NULL;
249: PetscInt d, coordSize;
250: PetscScalar *coords, *coordsIn = NULL;
251: PetscSection coordSection;
252: Vec coordinates;
254: PetscFunctionBegin;
255: PetscCallMPI(MPI_Comm_rank(comm, &rank));
257: if (rank == 0) {
258: FluentSection s;
259: do {
260: PetscCall(DMPlexCreateFluent_ReadSection(viewer, &s));
261: if (s.index == 2) { /* Dimension */
262: dim = s.nd;
264: } else if (s.index == 10 || s.index == 2010) { /* Vertices */
265: if (s.zoneID == 0) numVertices = s.last;
266: else {
267: PetscCheck(!coordsIn, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files");
268: coordsIn = (PetscScalar *)s.data;
269: }
271: } else if (s.index == 12 || s.index == 2012) { /* Cells */
272: if (s.zoneID == 0) numCells = s.last;
273: else {
274: switch (s.nd) {
275: case 0:
276: numCellVertices = PETSC_DETERMINE;
277: break;
278: case 1:
279: numCellVertices = 3;
280: break; /* triangular */
281: case 2:
282: numCellVertices = 4;
283: break; /* tetrahedral */
284: case 3:
285: numCellVertices = 4;
286: break; /* quadrilateral */
287: case 4:
288: numCellVertices = 8;
289: break; /* hexahedral */
290: case 5:
291: numCellVertices = 5;
292: break; /* pyramid */
293: case 6:
294: numCellVertices = 6;
295: break; /* wedge */
296: default:
297: numCellVertices = PETSC_DETERMINE;
298: }
299: }
301: } else if (s.index == 13 || s.index == 2013) { /* Facets */
302: if (s.zoneID == 0) { /* Header section */
303: numFaces = (PetscInt)(s.last - s.first + 1);
304: if (s.nd == 0 || s.nd == 5) numFaceVertices = PETSC_DETERMINE;
305: else numFaceVertices = s.nd;
306: } else { /* Data section */
307: unsigned int z;
309: PetscCheck(numFaceVertices == PETSC_DETERMINE || s.nd == numFaceVertices, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed facets in Fluent files are not supported");
310: PetscCheck(numFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file");
311: if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd;
312: numFaceEntries = numFaceVertices + 2;
313: if (!faces) PetscCall(PetscMalloc1(numFaces * numFaceEntries, &faces));
314: if (!faceZoneIDs) PetscCall(PetscMalloc1(numFaces, &faceZoneIDs));
315: PetscCall(PetscMemcpy(&faces[(s.first - 1) * numFaceEntries], s.data, (s.last - s.first + 1) * numFaceEntries * sizeof(PetscInt)));
316: /* Record the zoneID for each face set */
317: for (z = s.first - 1; z < s.last; z++) faceZoneIDs[z] = s.zoneID;
318: PetscCall(PetscFree(s.data));
319: }
320: }
321: } while (s.index >= 0);
322: }
323: PetscCallMPI(MPI_Bcast(&dim, 1, MPIU_INT, 0, comm));
324: PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension");
326: /* Allocate cell-vertex mesh */
327: PetscCall(DMCreate(comm, dm));
328: PetscCall(DMSetType(*dm, DMPLEX));
329: PetscCall(DMSetDimension(*dm, dim));
330: PetscCall(DMPlexSetChart(*dm, 0, numCells + numVertices));
331: if (rank == 0) {
332: PetscCheck(numCells >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown number of cells in Fluent file");
333: /* If no cell type was given we assume simplices */
334: if (numCellVertices == PETSC_DETERMINE) numCellVertices = numFaceVertices + 1;
335: for (c = 0; c < numCells; ++c) PetscCall(DMPlexSetConeSize(*dm, c, numCellVertices));
336: }
337: PetscCall(DMSetUp(*dm));
339: if (rank == 0 && faces) {
340: /* Derive cell-vertex list from face-vertex and face-cell maps */
341: PetscCall(PetscMalloc1(numCells * numCellVertices, &cellVertices));
342: for (c = 0; c < numCells * numCellVertices; c++) cellVertices[c] = -1;
343: for (f = 0; f < numFaces; f++) {
344: PetscInt *cell;
345: const PetscInt cl = faces[f * numFaceEntries + numFaceVertices];
346: const PetscInt cr = faces[f * numFaceEntries + numFaceVertices + 1];
347: const PetscInt *face = &(faces[f * numFaceEntries]);
349: if (cl > 0) {
350: cell = &(cellVertices[(cl - 1) * numCellVertices]);
351: for (v = 0; v < numFaceVertices; v++) {
352: PetscBool found = PETSC_FALSE;
353: for (c = 0; c < numCellVertices; c++) {
354: if (cell[c] < 0) break;
355: if (cell[c] == face[v] - 1 + numCells) {
356: found = PETSC_TRUE;
357: break;
358: }
359: }
360: if (!found) cell[c] = face[v] - 1 + numCells;
361: }
362: }
363: if (cr > 0) {
364: cell = &(cellVertices[(cr - 1) * numCellVertices]);
365: for (v = 0; v < numFaceVertices; v++) {
366: PetscBool found = PETSC_FALSE;
367: for (c = 0; c < numCellVertices; c++) {
368: if (cell[c] < 0) break;
369: if (cell[c] == face[v] - 1 + numCells) {
370: found = PETSC_TRUE;
371: break;
372: }
373: }
374: if (!found) cell[c] = face[v] - 1 + numCells;
375: }
376: }
377: }
378: for (c = 0; c < numCells; c++) PetscCall(DMPlexSetCone(*dm, c, &(cellVertices[c * numCellVertices])));
379: }
380: PetscCall(DMPlexSymmetrize(*dm));
381: PetscCall(DMPlexStratify(*dm));
382: if (interpolate) {
383: DM idm;
385: PetscCall(DMPlexInterpolate(*dm, &idm));
386: PetscCall(DMDestroy(dm));
387: *dm = idm;
388: }
390: if (rank == 0 && faces) {
391: PetscInt fi, joinSize, meetSize, *fverts, cells[2];
392: const PetscInt *join, *meet;
393: PetscCall(PetscMalloc1(numFaceVertices, &fverts));
394: /* Mark facets by finding the full join of all adjacent vertices */
395: for (f = 0; f < numFaces; f++) {
396: const PetscInt cl = faces[f * numFaceEntries + numFaceVertices] - 1;
397: const PetscInt cr = faces[f * numFaceEntries + numFaceVertices + 1] - 1;
398: if (cl > 0 && cr > 0) {
399: /* If we know both adjoining cells we can use a single-level meet */
400: cells[0] = cl;
401: cells[1] = cr;
402: PetscCall(DMPlexGetMeet(*dm, 2, cells, &meetSize, &meet));
403: PetscCheck(meetSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %" PetscInt_FMT, f);
404: PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", meet[0], faceZoneIDs[f]));
405: PetscCall(DMPlexRestoreMeet(*dm, numFaceVertices, fverts, &meetSize, &meet));
406: } else {
407: for (fi = 0; fi < numFaceVertices; fi++) fverts[fi] = faces[f * numFaceEntries + fi] + numCells - 1;
408: PetscCall(DMPlexGetFullJoin(*dm, numFaceVertices, fverts, &joinSize, &join));
409: PetscCheck(joinSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %" PetscInt_FMT, f);
410: PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], faceZoneIDs[f]));
411: PetscCall(DMPlexRestoreJoin(*dm, numFaceVertices, fverts, &joinSize, &join));
412: }
413: }
414: PetscCall(PetscFree(fverts));
415: }
417: { /* Create Face Sets label at all processes */
418: enum {
419: n = 1
420: };
421: PetscBool flag[n];
423: flag[0] = faceSets ? PETSC_TRUE : PETSC_FALSE;
424: PetscCallMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm));
425: if (flag[0]) PetscCall(DMCreateLabel(*dm, "Face Sets"));
426: }
428: /* Read coordinates */
429: PetscCall(DMGetCoordinateSection(*dm, &coordSection));
430: PetscCall(PetscSectionSetNumFields(coordSection, 1));
431: PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dim));
432: PetscCall(PetscSectionSetChart(coordSection, numCells, numCells + numVertices));
433: for (v = numCells; v < numCells + numVertices; ++v) {
434: PetscCall(PetscSectionSetDof(coordSection, v, dim));
435: PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dim));
436: }
437: PetscCall(PetscSectionSetUp(coordSection));
438: PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize));
439: PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates));
440: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
441: PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE));
442: PetscCall(VecSetType(coordinates, VECSTANDARD));
443: PetscCall(VecGetArray(coordinates, &coords));
444: if (rank == 0 && coordsIn) {
445: for (v = 0; v < numVertices; ++v) {
446: for (d = 0; d < dim; ++d) coords[v * dim + d] = coordsIn[v * dim + d];
447: }
448: }
449: PetscCall(VecRestoreArray(coordinates, &coords));
450: PetscCall(DMSetCoordinatesLocal(*dm, coordinates));
451: PetscCall(VecDestroy(&coordinates));
453: if (rank == 0) {
454: PetscCall(PetscFree(cellVertices));
455: PetscCall(PetscFree(faces));
456: PetscCall(PetscFree(faceZoneIDs));
457: PetscCall(PetscFree(coordsIn));
458: }
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }