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