Actual source code: plexfluent.c
petsc-3.9.4 2018-09-11
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;
22: PetscErrorCode ierr;
25: /* Create file viewer and build plex */
26: PetscViewerCreate(comm, &viewer);
27: PetscViewerSetType(viewer, PETSCVIEWERASCII);
28: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
29: PetscViewerFileSetName(viewer, filename);
30: DMPlexCreateFluent(comm, viewer, interpolate, dm);
31: PetscViewerDestroy(&viewer);
32: return(0);
33: }
35: static PetscErrorCode DMPlexCreateFluent_ReadString(PetscViewer viewer, char *buffer, char delim)
36: {
37: PetscInt ret, i = 0;
41: do {PetscViewerRead(viewer, &(buffer[i++]), 1, &ret, PETSC_CHAR);}
42: while (ret > 0 && buffer[i-1] != '\0' && buffer[i-1] != delim);
43: 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;
55: if (binary) {
56: /* Extract raw file descriptor to read binary block */
57: PetscViewerASCIIGetPointer(viewer, &file);
58: fflush(file); fdes = fileno(file);
59: }
61: if (!binary && dtype == PETSC_INT) {
62: char cbuf[256];
63: unsigned int ibuf;
64: int snum;
65: /* Parse hexadecimal ascii integers */
66: for (i = 0; i < count; i++) {
67: PetscViewerRead(viewer, cbuf, 1, NULL, PETSC_STRING);
68: snum = sscanf(cbuf, "%x", &ibuf);
69: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
70: ((PetscInt*)data)[i] = (PetscInt)ibuf;
71: }
72: } else if (binary && dtype == PETSC_INT) {
73: /* Always read 32-bit ints and cast to PetscInt */
74: int *ibuf;
75: PetscMalloc1(count, &ibuf);
76: PetscBinaryRead(fdes, ibuf, count, PETSC_ENUM);
77: PetscByteSwap(ibuf, PETSC_ENUM, count);
78: for (i = 0; i < count; i++) ((PetscInt*)data)[i] = (PetscInt)(ibuf[i]);
79: PetscFree(ibuf);
81: } else if (binary && dtype == PETSC_SCALAR) {
82: float *fbuf;
83: /* Always read 32-bit floats and cast to PetscScalar */
84: PetscMalloc1(count, &fbuf);
85: PetscBinaryRead(fdes, fbuf, count, PETSC_FLOAT);
86: PetscByteSwap(fbuf, PETSC_FLOAT, count);
87: for (i = 0; i < count; i++) ((PetscScalar*)data)[i] = (PetscScalar)(fbuf[i]);
88: PetscFree(fbuf);
89: } else {
90: PetscViewerASCIIRead(viewer, data, count, NULL, dtype);
91: }
92: return(0);
93: }
95: static PetscErrorCode DMPlexCreateFluent_ReadSection(PetscViewer viewer, FluentSection *s)
96: {
97: char buffer[PETSC_MAX_PATH_LEN];
98: int snum;
102: /* Fast-forward to next section and derive its index */
103: DMPlexCreateFluent_ReadString(viewer, buffer, '(');
104: DMPlexCreateFluent_ReadString(viewer, buffer, ' ');
105: snum = sscanf(buffer, "%d", &(s->index));
106: /* If we can't match an index return -1 to signal end-of-file */
107: if (snum < 1) {s->index = -1; return(0);}
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));
115: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
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));
120: if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
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));
133: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
134: if (s->zoneID == 0) { /* Header section */
135: snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
136: if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
137: } else { /* Data section */
138: snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
139: if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
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));
155: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
156: if (s->zoneID == 0) { /* Header section */
157: snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
158: if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
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));
162: if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
163: DMPlexCreateFluent_ReadString(viewer, buffer, '(');
164: switch (s->nd) {
165: case 0: numEntries = PETSC_DETERMINE; break;
166: case 2: numEntries = 2 + 2; break; /* linear */
167: case 3: numEntries = 2 + 3; break; /* triangular */
168: case 4: numEntries = 2 + 4; break; /* quadrilateral */
169: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file");
170: }
171: numFaces = s->last-s->first + 1;
172: if (numEntries != PETSC_DETERMINE) {
173: /* Allocate space only if we already know the size of the block */
174: PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);
175: }
176: for (f = 0; f < numFaces; f++) {
177: if (s->nd == 0) {
178: /* Determine the size of the block for "mixed" facets */
179: PetscInt numFaceVert = 0;
180: DMPlexCreateFluent_ReadValues(viewer, &numFaceVert, 1, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);
181: if (numEntries == PETSC_DETERMINE) {
182: numEntries = numFaceVert + 2;
183: PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);
184: } else {
185: if (numEntries != numFaceVert + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for mixed faces in Fluent files");
186: }
187: }
188: DMPlexCreateFluent_ReadValues(viewer, &(((PetscInt*)s->data)[f*numEntries]), numEntries, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);
189: }
190: s->nd = numEntries - 2;
191: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
192: }
193: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
195: } else { /* Unknown section type */
196: PetscInt depth = 1;
197: do {
198: /* Match parentheses when parsing unknown sections */
199: do {PetscViewerRead(viewer, &(buffer[0]), 1, NULL, PETSC_CHAR);}
200: while (buffer[0] != '(' && buffer[0] != ')');
201: if (buffer[0] == '(') depth++;
202: if (buffer[0] == ')') depth--;
203: } while (depth > 0);
204: DMPlexCreateFluent_ReadString(viewer, buffer, '\n');
205: }
206: return(0);
207: }
209: /*@C
210: DMPlexCreateFluent - Create a DMPlex mesh from a Fluent mesh file.
212: Collective on comm
214: Input Parameters:
215: + comm - The MPI communicator
216: . viewer - The Viewer associated with a Fluent mesh file
217: - interpolate - Create faces and edges in the mesh
219: Output Parameter:
220: . dm - The DM object representing the mesh
222: Note: http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm
224: Level: beginner
226: .keywords: mesh, fluent, case
227: .seealso: DMPLEX, DMCreate()
228: @*/
229: PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
230: {
231: PetscMPIInt rank;
232: PetscInt c, v, dim = PETSC_DETERMINE, numCells = 0, numVertices = 0, numCellVertices = PETSC_DETERMINE;
233: PetscInt numFaces = PETSC_DETERMINE, f, numFaceEntries = PETSC_DETERMINE, numFaceVertices = PETSC_DETERMINE;
234: PetscInt *faces = NULL, *cellVertices = NULL, *faceZoneIDs = NULL;
235: PetscInt d, coordSize;
236: PetscScalar *coords, *coordsIn = NULL;
237: PetscSection coordSection;
238: Vec coordinates;
242: MPI_Comm_rank(comm, &rank);
244: if (!rank) {
245: FluentSection s;
246: do {
247: DMPlexCreateFluent_ReadSection(viewer, &s);
248: if (s.index == 2) { /* Dimension */
249: dim = s.nd;
251: } else if (s.index == 10 || s.index == 2010) { /* Vertices */
252: if (s.zoneID == 0) numVertices = s.last;
253: else {
254: if (coordsIn) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files");
255: coordsIn = (PetscScalar *) s.data;
256: }
258: } else if (s.index == 12 || s.index == 2012) { /* Cells */
259: if (s.zoneID == 0) numCells = s.last;
260: else {
261: switch (s.nd) {
262: case 0: numCellVertices = PETSC_DETERMINE; break;
263: case 1: numCellVertices = 3; break; /* triangular */
264: case 2: numCellVertices = 4; break; /* tetrahedral */
265: case 3: numCellVertices = 4; break; /* quadrilateral */
266: case 4: numCellVertices = 8; break; /* hexahedral */
267: case 5: numCellVertices = 5; break; /* pyramid */
268: case 6: numCellVertices = 6; break; /* wedge */
269: default: numCellVertices = PETSC_DETERMINE;
270: }
271: }
273: } else if (s.index == 13 || s.index == 2013) { /* Facets */
274: if (s.zoneID == 0) { /* Header section */
275: numFaces = (PetscInt) (s.last - s.first + 1);
276: if (s.nd == 0 || s.nd == 5) numFaceVertices = PETSC_DETERMINE;
277: else numFaceVertices = s.nd;
278: } else { /* Data section */
279: unsigned int z;
281: if (numFaceVertices != PETSC_DETERMINE && s.nd != numFaceVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed facets in Fluent files are not supported");
282: if (numFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file");
283: if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd;
284: numFaceEntries = numFaceVertices + 2;
285: if (!faces) {PetscMalloc1(numFaces*numFaceEntries, &faces);}
286: if (!faceZoneIDs) {PetscMalloc1(numFaces, &faceZoneIDs);}
287: PetscMemcpy(&(faces[(s.first-1)*numFaceEntries]), s.data, (s.last-s.first+1)*numFaceEntries*sizeof(PetscInt));
288: /* Record the zoneID for each face set */
289: for (z = s.first -1; z < s.last; z++) faceZoneIDs[z] = s.zoneID;
290: PetscFree(s.data);
291: }
292: }
293: } while (s.index >= 0);
294: }
295: MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
296: if (dim < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension");
298: /* Allocate cell-vertex mesh */
299: DMCreate(comm, dm);
300: DMSetType(*dm, DMPLEX);
301: DMSetDimension(*dm, dim);
302: DMPlexSetChart(*dm, 0, numCells + numVertices);
303: if (!rank) {
304: if (numCells < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown number of cells in Fluent file");
305: /* If no cell type was given we assume simplices */
306: if (numCellVertices == PETSC_DETERMINE) numCellVertices = numFaceVertices + 1;
307: for (c = 0; c < numCells; ++c) {DMPlexSetConeSize(*dm, c, numCellVertices);}
308: }
309: DMSetUp(*dm);
311: if (!rank && faces) {
312: /* Derive cell-vertex list from face-vertex and face-cell maps */
313: PetscMalloc1(numCells*numCellVertices, &cellVertices);
314: for (c = 0; c < numCells*numCellVertices; c++) cellVertices[c] = -1;
315: for (f = 0; f < numFaces; f++) {
316: PetscInt *cell;
317: const PetscInt cl = faces[f*numFaceEntries + numFaceVertices];
318: const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1];
319: const PetscInt *face = &(faces[f*numFaceEntries]);
321: if (cl > 0) {
322: cell = &(cellVertices[(cl-1) * numCellVertices]);
323: for (v = 0; v < numFaceVertices; v++) {
324: PetscBool found = PETSC_FALSE;
325: for (c = 0; c < numCellVertices; c++) {
326: if (cell[c] < 0) break;
327: if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
328: }
329: if (!found) cell[c] = face[v]-1 + numCells;
330: }
331: }
332: if (cr > 0) {
333: cell = &(cellVertices[(cr-1) * numCellVertices]);
334: for (v = 0; v < numFaceVertices; v++) {
335: PetscBool found = PETSC_FALSE;
336: for (c = 0; c < numCellVertices; c++) {
337: if (cell[c] < 0) break;
338: if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
339: }
340: if (!found) cell[c] = face[v]-1 + numCells;
341: }
342: }
343: }
344: for (c = 0; c < numCells; c++) {
345: DMPlexSetCone(*dm, c, &(cellVertices[c*numCellVertices]));
346: }
347: }
348: DMPlexSymmetrize(*dm);
349: DMPlexStratify(*dm);
350: if (interpolate) {
351: DM idm;
353: DMPlexInterpolate(*dm, &idm);
354: DMDestroy(dm);
355: *dm = idm;
356: }
358: if (!rank && faces) {
359: PetscInt fi, joinSize, meetSize, *fverts, cells[2];
360: const PetscInt *join, *meet;
361: PetscMalloc1(numFaceVertices, &fverts);
362: /* Mark facets by finding the full join of all adjacent vertices */
363: for (f = 0; f < numFaces; f++) {
364: const PetscInt cl = faces[f*numFaceEntries + numFaceVertices] - 1;
365: const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1] - 1;
366: if (cl > 0 && cr > 0) {
367: /* If we know both adjoining cells we can use a single-level meet */
368: cells[0] = cl; cells[1] = cr;
369: DMPlexGetMeet(*dm, 2, cells, &meetSize, &meet);
370: if (meetSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %d", f);
371: DMSetLabelValue(*dm, "Face Sets", meet[0], faceZoneIDs[f]);
372: DMPlexRestoreMeet(*dm, numFaceVertices, fverts, &meetSize, &meet);
373: } else {
374: for (fi = 0; fi < numFaceVertices; fi++) fverts[fi] = faces[f*numFaceEntries + fi] + numCells - 1;
375: DMPlexGetFullJoin(*dm, numFaceVertices, fverts, &joinSize, &join);
376: if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %d", f);
377: DMSetLabelValue(*dm, "Face Sets", join[0], faceZoneIDs[f]);
378: DMPlexRestoreJoin(*dm, numFaceVertices, fverts, &joinSize, &join);
379: }
380: }
381: PetscFree(fverts);
382: }
384: /* Read coordinates */
385: DMGetCoordinateSection(*dm, &coordSection);
386: PetscSectionSetNumFields(coordSection, 1);
387: PetscSectionSetFieldComponents(coordSection, 0, dim);
388: PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
389: for (v = numCells; v < numCells+numVertices; ++v) {
390: PetscSectionSetDof(coordSection, v, dim);
391: PetscSectionSetFieldDof(coordSection, v, 0, dim);
392: }
393: PetscSectionSetUp(coordSection);
394: PetscSectionGetStorageSize(coordSection, &coordSize);
395: VecCreate(PETSC_COMM_SELF, &coordinates);
396: PetscObjectSetName((PetscObject) coordinates, "coordinates");
397: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
398: VecSetType(coordinates, VECSTANDARD);
399: VecGetArray(coordinates, &coords);
400: if (!rank && coordsIn) {
401: for (v = 0; v < numVertices; ++v) {
402: for (d = 0; d < dim; ++d) {
403: coords[v*dim+d] = coordsIn[v*dim+d];
404: }
405: }
406: }
407: VecRestoreArray(coordinates, &coords);
408: DMSetCoordinatesLocal(*dm, coordinates);
409: VecDestroy(&coordinates);
410: if (!rank) {
411: if (cellVertices) {PetscFree(cellVertices);}
412: PetscFree(faces);
413: PetscFree(faceZoneIDs);
414: PetscFree(coordsIn);
415: }
416: return(0);
417: }