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