Actual source code: plexfluent.c
petsc-3.8.4 2018-03-24
1: #define PETSCDM_DLL
2: #include <petsc/private/dmpleximpl.h>
4: /*@C
5: DMPlexCreateFluentFromFile - Create a DMPlex mesh from a Fluent mesh file
7: + comm - The MPI communicator
8: . filename - Name of the Fluent mesh file
9: - interpolate - Create faces and edges in the mesh
11: Output Parameter:
12: . dm - The DM object representing the mesh
14: Level: beginner
16: .seealso: DMPlexCreateFromFile(), DMPlexCreateFluent(), DMPlexCreate()
17: @*/
18: PetscErrorCode DMPlexCreateFluentFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
19: {
20: PetscViewer viewer;
21: PetscErrorCode ierr;
24: /* Create file viewer and build plex */
25: PetscViewerCreate(comm, &viewer);
26: PetscViewerSetType(viewer, PETSCVIEWERASCII);
27: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
28: PetscViewerFileSetName(viewer, filename);
29: DMPlexCreateFluent(comm, viewer, interpolate, dm);
30: PetscViewerDestroy(&viewer);
31: return(0);
32: }
34: static PetscErrorCode DMPlexCreateFluent_ReadString(PetscViewer viewer, char *buffer, char delim)
35: {
36: 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: buffer[i] = '\0';
43: return(0);
44: }
46: static PetscErrorCode DMPlexCreateFluent_ReadValues(PetscViewer viewer, void *data, PetscInt count, PetscDataType dtype, PetscBool binary)
47: {
48: int fdes=0;
49: FILE *file;
50: PetscInt i;
54: if (binary) {
55: /* Extract raw file descriptor to read binary block */
56: PetscViewerASCIIGetPointer(viewer, &file);
57: fflush(file); 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);
68: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
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, 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, 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;
101: /* Fast-forward to next section and derive its index */
102: DMPlexCreateFluent_ReadString(viewer, buffer, '(');
103: DMPlexCreateFluent_ReadString(viewer, buffer, ' ');
104: snum = sscanf(buffer, "%d", &(s->index));
105: /* If we can't match an index return -1 to signal end-of-file */
106: if (snum < 1) {s->index = -1; return(0);}
108: if (s->index == 0) { /* Comment */
109: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
111: } else if (s->index == 2) { /* Dimension */
112: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
113: snum = sscanf(buffer, "%d", &(s->nd));
114: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
116: } else if (s->index == 10 || s->index == 2010) { /* Vertices */
117: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
118: snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
119: if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
120: if (s->zoneID > 0) {
121: PetscInt numCoords = s->last - s->first + 1;
122: DMPlexCreateFluent_ReadString(viewer, buffer, '(');
123: PetscMalloc1(s->nd*numCoords, (PetscScalar**)&s->data);
124: DMPlexCreateFluent_ReadValues(viewer, s->data, s->nd*numCoords, PETSC_SCALAR, s->index==2010 ? PETSC_TRUE : PETSC_FALSE);
125: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
126: }
127: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
129: } else if (s->index == 12 || s->index == 2012) { /* Cells */
130: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
131: snum = sscanf(buffer, "(%x", &(s->zoneID));
132: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
133: if (s->zoneID == 0) { /* Header section */
134: snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
135: if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
136: } else { /* Data section */
137: snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
138: if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
139: if (s->nd == 0) {
140: /* Read cell type definitions for mixed cells */
141: PetscInt numCells = s->last - s->first + 1;
142: DMPlexCreateFluent_ReadString(viewer, buffer, '(');
143: PetscMalloc1(numCells, (PetscInt**)&s->data);
144: DMPlexCreateFluent_ReadValues(viewer, s->data, numCells, PETSC_INT, s->index==2012 ? PETSC_TRUE : PETSC_FALSE);
145: PetscFree(s->data);
146: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
147: }
148: }
149: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
151: } else if (s->index == 13 || s->index == 2013) { /* Faces */
152: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
153: snum = sscanf(buffer, "(%x", &(s->zoneID));
154: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
155: if (s->zoneID == 0) { /* Header section */
156: snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
157: if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
158: } else { /* Data section */
159: PetscInt f, numEntries, numFaces;
160: snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
161: if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
162: DMPlexCreateFluent_ReadString(viewer, buffer, '(');
163: switch (s->nd) {
164: case 0: numEntries = PETSC_DETERMINE; break;
165: case 2: numEntries = 2 + 2; break; /* linear */
166: case 3: numEntries = 2 + 3; break; /* triangular */
167: case 4: numEntries = 2 + 4; break; /* quadrilateral */
168: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file");
169: }
170: numFaces = s->last-s->first + 1;
171: if (numEntries != PETSC_DETERMINE) {
172: /* Allocate space only if we already know the size of the block */
173: PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);
174: }
175: for (f = 0; f < numFaces; f++) {
176: if (s->nd == 0) {
177: /* Determine the size of the block for "mixed" facets */
178: PetscInt numFaceVert = 0;
179: DMPlexCreateFluent_ReadValues(viewer, &numFaceVert, 1, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);
180: if (numEntries == PETSC_DETERMINE) {
181: numEntries = numFaceVert + 2;
182: PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);
183: } else {
184: if (numEntries != numFaceVert + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for mixed faces in Fluent files");
185: }
186: }
187: DMPlexCreateFluent_ReadValues(viewer, &(((PetscInt*)s->data)[f*numEntries]), numEntries, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);
188: }
189: s->nd = numEntries - 2;
190: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
191: }
192: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
194: } else { /* Unknown section type */
195: PetscInt depth = 1;
196: do {
197: /* Match parentheses when parsing unknown sections */
198: do {PetscViewerRead(viewer, &(buffer[0]), 1, NULL, PETSC_CHAR);}
199: while (buffer[0] != '(' && buffer[0] != ')');
200: if (buffer[0] == '(') depth++;
201: if (buffer[0] == ')') depth--;
202: } while (depth > 0);
203: DMPlexCreateFluent_ReadString(viewer, buffer, '\n');
204: }
205: return(0);
206: }
208: /*@C
209: DMPlexCreateFluent - Create a DMPlex mesh from a Fluent mesh file.
211: Collective on comm
213: Input Parameters:
214: + comm - The MPI communicator
215: . viewer - The Viewer associated with a Fluent mesh file
216: - interpolate - Create faces and edges in the mesh
218: Output Parameter:
219: . dm - The DM object representing the mesh
221: Note: http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm
223: Level: beginner
225: .keywords: mesh, fluent, case
226: .seealso: DMPLEX, DMCreate()
227: @*/
228: PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
229: {
230: PetscMPIInt rank;
231: PetscInt c, v, dim = PETSC_DETERMINE, numCells = 0, numVertices = 0, numCellVertices = PETSC_DETERMINE;
232: PetscInt numFaces = PETSC_DETERMINE, f, numFaceEntries = PETSC_DETERMINE, numFaceVertices = PETSC_DETERMINE;
233: PetscInt *faces = NULL, *cellVertices = NULL, *faceZoneIDs = NULL;
234: PetscInt d, coordSize;
235: PetscScalar *coords, *coordsIn = NULL;
236: PetscSection coordSection;
237: Vec coordinates;
241: MPI_Comm_rank(comm, &rank);
243: if (!rank) {
244: FluentSection s;
245: do {
246: DMPlexCreateFluent_ReadSection(viewer, &s);
247: if (s.index == 2) { /* Dimension */
248: dim = s.nd;
250: } else if (s.index == 10 || s.index == 2010) { /* Vertices */
251: if (s.zoneID == 0) numVertices = s.last;
252: else {
253: if (coordsIn) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files");
254: coordsIn = (PetscScalar *) s.data;
255: }
257: } else if (s.index == 12 || s.index == 2012) { /* Cells */
258: if (s.zoneID == 0) numCells = s.last;
259: else {
260: switch (s.nd) {
261: case 0: numCellVertices = PETSC_DETERMINE; break;
262: case 1: numCellVertices = 3; break; /* triangular */
263: case 2: numCellVertices = 4; break; /* tetrahedral */
264: case 3: numCellVertices = 4; break; /* quadrilateral */
265: case 4: numCellVertices = 8; break; /* hexahedral */
266: case 5: numCellVertices = 5; break; /* pyramid */
267: case 6: numCellVertices = 6; break; /* wedge */
268: default: numCellVertices = PETSC_DETERMINE;
269: }
270: }
272: } else if (s.index == 13 || s.index == 2013) { /* Facets */
273: if (s.zoneID == 0) { /* Header section */
274: numFaces = (PetscInt) (s.last - s.first + 1);
275: if (s.nd == 0 || s.nd == 5) numFaceVertices = PETSC_DETERMINE;
276: else numFaceVertices = s.nd;
277: } else { /* Data section */
278: unsigned int z;
280: if (numFaceVertices != PETSC_DETERMINE && s.nd != numFaceVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed facets in Fluent files are not supported");
281: if (numFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file");
282: if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd;
283: numFaceEntries = numFaceVertices + 2;
284: if (!faces) {PetscMalloc1(numFaces*numFaceEntries, &faces);}
285: if (!faceZoneIDs) {PetscMalloc1(numFaces, &faceZoneIDs);}
286: PetscMemcpy(&(faces[(s.first-1)*numFaceEntries]), s.data, (s.last-s.first+1)*numFaceEntries*sizeof(PetscInt));
287: /* Record the zoneID for each face set */
288: for (z = s.first -1; z < s.last; z++) faceZoneIDs[z] = s.zoneID;
289: PetscFree(s.data);
290: }
291: }
292: } while (s.index >= 0);
293: }
294: MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
295: if (dim < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension");
297: /* Allocate cell-vertex mesh */
298: DMCreate(comm, dm);
299: DMSetType(*dm, DMPLEX);
300: DMSetDimension(*dm, dim);
301: DMPlexSetChart(*dm, 0, numCells + numVertices);
302: if (!rank) {
303: if (numCells < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown number of cells in Fluent file");
304: /* If no cell type was given we assume simplices */
305: if (numCellVertices == PETSC_DETERMINE) numCellVertices = numFaceVertices + 1;
306: for (c = 0; c < numCells; ++c) {DMPlexSetConeSize(*dm, c, numCellVertices);}
307: }
308: DMSetUp(*dm);
310: if (!rank && faces) {
311: /* Derive cell-vertex list from face-vertex and face-cell maps */
312: PetscMalloc1(numCells*numCellVertices, &cellVertices);
313: for (c = 0; c < numCells*numCellVertices; c++) cellVertices[c] = -1;
314: for (f = 0; f < numFaces; f++) {
315: PetscInt *cell;
316: const PetscInt cl = faces[f*numFaceEntries + numFaceVertices];
317: const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1];
318: const PetscInt *face = &(faces[f*numFaceEntries]);
320: if (cl > 0) {
321: cell = &(cellVertices[(cl-1) * numCellVertices]);
322: for (v = 0; v < numFaceVertices; v++) {
323: PetscBool found = PETSC_FALSE;
324: for (c = 0; c < numCellVertices; c++) {
325: if (cell[c] < 0) break;
326: if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
327: }
328: if (!found) cell[c] = face[v]-1 + numCells;
329: }
330: }
331: if (cr > 0) {
332: cell = &(cellVertices[(cr-1) * numCellVertices]);
333: for (v = 0; v < numFaceVertices; v++) {
334: PetscBool found = PETSC_FALSE;
335: for (c = 0; c < numCellVertices; c++) {
336: if (cell[c] < 0) break;
337: if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
338: }
339: if (!found) cell[c] = face[v]-1 + numCells;
340: }
341: }
342: }
343: for (c = 0; c < numCells; c++) {
344: DMPlexSetCone(*dm, c, &(cellVertices[c*numCellVertices]));
345: }
346: }
347: DMPlexSymmetrize(*dm);
348: DMPlexStratify(*dm);
349: if (interpolate) {
350: DM idm = NULL;
352: DMPlexInterpolate(*dm, &idm);
353: DMDestroy(dm);
354: *dm = idm;
355: }
357: if (!rank && faces) {
358: PetscInt fi, joinSize, meetSize, *fverts, cells[2];
359: const PetscInt *join, *meet;
360: PetscMalloc1(numFaceVertices, &fverts);
361: /* Mark facets by finding the full join of all adjacent vertices */
362: for (f = 0; f < numFaces; f++) {
363: const PetscInt cl = faces[f*numFaceEntries + numFaceVertices] - 1;
364: const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1] - 1;
365: if (cl > 0 && cr > 0) {
366: /* If we know both adjoining cells we can use a single-level meet */
367: cells[0] = cl; cells[1] = cr;
368: DMPlexGetMeet(*dm, 2, cells, &meetSize, &meet);
369: if (meetSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %d", f);
370: DMSetLabelValue(*dm, "Face Sets", meet[0], faceZoneIDs[f]);
371: DMPlexRestoreMeet(*dm, numFaceVertices, fverts, &meetSize, &meet);
372: } else {
373: for (fi = 0; fi < numFaceVertices; fi++) fverts[fi] = faces[f*numFaceEntries + fi] + numCells - 1;
374: DMPlexGetFullJoin(*dm, numFaceVertices, fverts, &joinSize, &join);
375: if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %d", f);
376: DMSetLabelValue(*dm, "Face Sets", join[0], faceZoneIDs[f]);
377: DMPlexRestoreJoin(*dm, numFaceVertices, fverts, &joinSize, &join);
378: }
379: }
380: PetscFree(fverts);
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 && 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);
409: if (!rank) {
410: if (cellVertices) {PetscFree(cellVertices);}
411: PetscFree(faces);
412: PetscFree(faceZoneIDs);
413: PetscFree(coordsIn);
414: }
415: return(0);
416: }