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