Actual source code: plexfluent.c
petsc-3.7.5 2017-01-01
1: #define PETSCDM_DLL
2: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
6: /*@C
7: DMPlexCreateFluentFromFile - Create a DMPlex mesh from a Fluent mesh file
9: + comm - The MPI communicator
10: . filename - Name of the Fluent mesh file
11: - interpolate - Create faces and edges in the mesh
13: Output Parameter:
14: . dm - The DM object representing the mesh
16: Level: beginner
18: .seealso: DMPlexCreateFromFile(), DMPlexCreateFluent(), DMPlexCreate()
19: @*/
20: PetscErrorCode DMPlexCreateFluentFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
21: {
22: PetscViewer viewer;
23: PetscErrorCode ierr;
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: }
38: PetscErrorCode DMPlexCreateFluent_ReadString(PetscViewer viewer, char *buffer, char delim)
39: {
40: PetscInt ret, i = 0;
44: do {PetscViewerRead(viewer, &(buffer[i++]), 1, &ret, PETSC_CHAR);}
45: while (ret > 0 && buffer[i-1] != '\0' && buffer[i-1] != delim);
46: buffer[i] = '\0';
47: return(0);
48: }
52: PetscErrorCode DMPlexCreateFluent_ReadValues(PetscViewer viewer, void *data, PetscInt count, PetscDataType dtype, PetscBool binary)
53: {
54: int fdes=0;
55: FILE *file;
56: PetscInt i;
60: if (binary) {
61: /* Extract raw file descriptor to read binary block */
62: PetscViewerASCIIGetPointer(viewer, &file);
63: fflush(file); fdes = fileno(file);
64: }
66: if (!binary && dtype == PETSC_INT) {
67: char cbuf[256];
68: int ibuf, snum;
69: /* Parse hexadecimal ascii integers */
70: for (i = 0; i < count; i++) {
71: PetscViewerRead(viewer, cbuf, 1, NULL, PETSC_STRING);
72: snum = sscanf(cbuf, "%x", &ibuf);
73: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
74: ((PetscInt*)data)[i] = (PetscInt)ibuf;
75: }
76: } else if (binary && dtype == PETSC_INT) {
77: /* Always read 32-bit ints and cast to PetscInt */
78: int *ibuf;
79: PetscMalloc1(count, &ibuf);
80: PetscBinaryRead(fdes, ibuf, count, PETSC_ENUM);
81: PetscByteSwap(ibuf, PETSC_ENUM, count);
82: for (i = 0; i < count; i++) ((PetscInt*)data)[i] = (PetscInt)(ibuf[i]);
83: PetscFree(ibuf);
85: } else if (binary && dtype == PETSC_SCALAR) {
86: float *fbuf;
87: /* Always read 32-bit floats and cast to PetscScalar */
88: PetscMalloc1(count, &fbuf);
89: PetscBinaryRead(fdes, fbuf, count, PETSC_FLOAT);
90: PetscByteSwap(fbuf, PETSC_FLOAT, count);
91: for (i = 0; i < count; i++) ((PetscScalar*)data)[i] = (PetscScalar)(fbuf[i]);
92: PetscFree(fbuf);
93: } else {
94: PetscViewerASCIIRead(viewer, data, count, NULL, dtype);
95: }
96: return(0);
97: }
101: PetscErrorCode DMPlexCreateFluent_ReadSection(PetscViewer viewer, FluentSection *s)
102: {
103: char buffer[PETSC_MAX_PATH_LEN];
104: int snum;
108: /* Fast-forward to next section and derive its index */
109: DMPlexCreateFluent_ReadString(viewer, buffer, '(');
110: DMPlexCreateFluent_ReadString(viewer, buffer, ' ');
111: snum = sscanf(buffer, "%d", &(s->index));
112: /* If we can't match an index return -1 to signal end-of-file */
113: if (snum < 1) {s->index = -1; return(0);}
115: if (s->index == 0) { /* Comment */
116: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
118: } else if (s->index == 2) { /* Dimension */
119: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
120: snum = sscanf(buffer, "%d", &(s->nd));
121: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
123: } else if (s->index == 10 || s->index == 2010) { /* Vertices */
124: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
125: snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
126: if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
127: if (s->zoneID > 0) {
128: PetscInt numCoords = s->last - s->first + 1;
129: DMPlexCreateFluent_ReadString(viewer, buffer, '(');
130: PetscMalloc1(s->nd*numCoords, (PetscScalar**)&s->data);
131: DMPlexCreateFluent_ReadValues(viewer, s->data, s->nd*numCoords, PETSC_SCALAR, s->index==2010 ? PETSC_TRUE : PETSC_FALSE);
132: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
133: }
134: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
136: } else if (s->index == 12 || s->index == 2012) { /* Cells */
137: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
138: snum = sscanf(buffer, "(%x", &(s->zoneID));
139: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
140: if (s->zoneID == 0) { /* Header section */
141: snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
142: if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
143: } else { /* Data section */
144: snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
145: if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
146: if (s->nd == 0) {
147: /* Read cell type definitions for mixed cells */
148: PetscInt numCells = s->last - s->first + 1;
149: DMPlexCreateFluent_ReadString(viewer, buffer, '(');
150: PetscMalloc1(numCells, (PetscInt**)&s->data);
151: DMPlexCreateFluent_ReadValues(viewer, s->data, numCells, PETSC_INT, s->index==2012 ? PETSC_TRUE : PETSC_FALSE);
152: PetscFree(s->data);
153: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
154: }
155: }
156: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
158: } else if (s->index == 13 || s->index == 2013) { /* Faces */
159: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
160: snum = sscanf(buffer, "(%x", &(s->zoneID));
161: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
162: if (s->zoneID == 0) { /* Header section */
163: snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
164: if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
165: } else { /* Data section */
166: PetscInt f, numEntries, numFaces, numFaceVert;
167: snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
168: if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
169: DMPlexCreateFluent_ReadString(viewer, buffer, '(');
170: switch (s->nd) {
171: case 0: numEntries = PETSC_DETERMINE; break;
172: case 2: numEntries = 2 + 2; break; /* linear */
173: case 3: numEntries = 2 + 3; break; /* triangular */
174: case 4: numEntries = 2 + 4; break; /* quadrilateral */
175: default: 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: DMPlexCreateFluent_ReadValues(viewer, &numFaceVert, 1, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);
186: if (numEntries == PETSC_DETERMINE) {
187: numEntries = numFaceVert + 2;
188: PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);
189: } else {
190: if (numEntries != numFaceVert + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for mixed faces in Fluent files");
191: }
192: }
193: DMPlexCreateFluent_ReadValues(viewer, &(((PetscInt*)s->data)[f*numEntries]), numEntries, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);
194: }
195: s->nd = numEntries - 2;
196: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
197: }
198: DMPlexCreateFluent_ReadString(viewer, buffer, ')');
200: } else { /* Unknown section type */
201: PetscInt depth = 1;
202: do {
203: /* Match parentheses when parsing unknown sections */
204: do {PetscViewerRead(viewer, &(buffer[0]), 1, NULL, PETSC_CHAR);}
205: while (buffer[0] != '(' && buffer[0] != ')');
206: if (buffer[0] == '(') depth++;
207: if (buffer[0] == ')') depth--;
208: } while (depth > 0);
209: DMPlexCreateFluent_ReadString(viewer, buffer, '\n');
210: }
211: return(0);
212: }
216: /*@C
217: DMPlexCreateFluent - Create a DMPlex mesh from a Fluent mesh file.
219: Collective on comm
221: Input Parameters:
222: + comm - The MPI communicator
223: . viewer - The Viewer associated with a Fluent mesh file
224: - interpolate - Create faces and edges in the mesh
226: Output Parameter:
227: . dm - The DM object representing the mesh
229: Note: http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm
231: Level: beginner
233: .keywords: mesh, fluent, case
234: .seealso: DMPLEX, DMCreate()
235: @*/
236: PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
237: {
238: PetscMPIInt rank;
239: PetscInt c, f, v, dim = PETSC_DETERMINE, numCells = 0, numVertices = 0, numCellVertices = PETSC_DETERMINE;
240: PetscInt numFaces = PETSC_DETERMINE, numFaceEntries = PETSC_DETERMINE, numFaceVertices = PETSC_DETERMINE;
241: PetscInt *faces = NULL, *cellVertices = NULL, *faceZoneIDs = NULL;
242: PetscInt d, coordSize;
243: PetscScalar *coords, *coordsIn = NULL;
244: PetscSection coordSection;
245: Vec coordinates;
249: MPI_Comm_rank(comm, &rank);
251: if (!rank) {
252: FluentSection s;
253: do {
254: DMPlexCreateFluent_ReadSection(viewer, &s);
255: if (s.index == 2) { /* Dimension */
256: dim = s.nd;
258: } else if (s.index == 10 || s.index == 2010) { /* Vertices */
259: if (s.zoneID == 0) numVertices = s.last;
260: else {
261: if (coordsIn) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files");
262: coordsIn = (PetscScalar *) s.data;
263: }
265: } else if (s.index == 12 || s.index == 2012) { /* Cells */
266: if (s.zoneID == 0) numCells = s.last;
267: else {
268: switch (s.nd) {
269: case 0: numCellVertices = PETSC_DETERMINE; break;
270: case 1: numCellVertices = 3; break; /* triangular */
271: case 2: numCellVertices = 4; break; /* tetrahedral */
272: case 3: numCellVertices = 4; break; /* quadrilateral */
273: case 4: numCellVertices = 8; break; /* hexahedral */
274: case 5: numCellVertices = 5; break; /* pyramid */
275: case 6: numCellVertices = 6; break; /* wedge */
276: default: numCellVertices = PETSC_DETERMINE;
277: }
278: }
280: } else if (s.index == 13 || s.index == 2013) { /* Facets */
281: if (s.zoneID == 0) { /* Header section */
282: numFaces = s.last - s.first + 1;
283: if (s.nd == 0 || s.nd == 5) numFaceVertices = PETSC_DETERMINE;
284: else numFaceVertices = s.nd;
285: } else { /* Data section */
286: if (numFaceVertices != PETSC_DETERMINE && s.nd != numFaceVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed facets in Fluent files are not supported");
287: if (numFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file");
288: if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd;
289: numFaceEntries = numFaceVertices + 2;
290: if (!faces) {PetscMalloc1(numFaces*numFaceEntries, &faces);}
291: if (!faceZoneIDs) {PetscMalloc1(numFaces, &faceZoneIDs);}
292: PetscMemcpy(&(faces[(s.first-1)*numFaceEntries]), s.data, (s.last-s.first+1)*numFaceEntries*sizeof(PetscInt));
293: /* Record the zoneID for each face set */
294: for (f = s.first -1; f < s.last; f++) faceZoneIDs[f] = s.zoneID;
295: PetscFree(s.data);
296: }
297: }
298: } while (s.index >= 0);
299: }
300: MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);
301: if (dim < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension");
303: /* Allocate cell-vertex mesh */
304: DMCreate(comm, dm);
305: DMSetType(*dm, DMPLEX);
306: DMSetDimension(*dm, dim);
307: DMPlexSetChart(*dm, 0, numCells + numVertices);
308: if (!rank) {
309: if (numCells < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown number of cells in Fluent file");
310: /* If no cell type was given we assume simplices */
311: if (numCellVertices == PETSC_DETERMINE) numCellVertices = numFaceVertices + 1;
312: for (c = 0; c < numCells; ++c) {DMPlexSetConeSize(*dm, c, numCellVertices);}
313: }
314: DMSetUp(*dm);
316: if (!rank && faces) {
317: /* Derive cell-vertex list from face-vertex and face-cell maps */
318: PetscMalloc1(numCells*numCellVertices, &cellVertices);
319: for (c = 0; c < numCells*numCellVertices; c++) cellVertices[c] = -1;
320: for (f = 0; f < numFaces; f++) {
321: PetscInt *cell;
322: const PetscInt cl = faces[f*numFaceEntries + numFaceVertices];
323: const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1];
324: const PetscInt *face = &(faces[f*numFaceEntries]);
326: if (cl > 0) {
327: cell = &(cellVertices[(cl-1) * numCellVertices]);
328: for (v = 0; v < numFaceVertices; v++) {
329: PetscBool found = PETSC_FALSE;
330: for (c = 0; c < numCellVertices; c++) {
331: if (cell[c] < 0) break;
332: if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
333: }
334: if (!found) cell[c] = face[v]-1 + numCells;
335: }
336: }
337: if (cr > 0) {
338: cell = &(cellVertices[(cr-1) * numCellVertices]);
339: for (v = 0; v < numFaceVertices; v++) {
340: PetscBool found = PETSC_FALSE;
341: for (c = 0; c < numCellVertices; c++) {
342: if (cell[c] < 0) break;
343: if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
344: }
345: if (!found) cell[c] = face[v]-1 + numCells;
346: }
347: }
348: }
349: for (c = 0; c < numCells; c++) {
350: DMPlexSetCone(*dm, c, &(cellVertices[c*numCellVertices]));
351: }
352: }
353: DMPlexSymmetrize(*dm);
354: DMPlexStratify(*dm);
355: if (interpolate) {
356: DM idm = NULL;
358: DMPlexInterpolate(*dm, &idm);
359: DMDestroy(dm);
360: *dm = idm;
361: }
363: if (!rank && faces) {
364: PetscInt fi, joinSize, meetSize, *fverts, cells[2];
365: const PetscInt *join, *meet;
366: PetscMalloc1(numFaceVertices, &fverts);
367: /* Mark facets by finding the full join of all adjacent vertices */
368: for (f = 0; f < numFaces; f++) {
369: const PetscInt cl = faces[f*numFaceEntries + numFaceVertices] - 1;
370: const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1] - 1;
371: if (cl > 0 && cr > 0) {
372: /* If we know both adjoining cells we can use a single-level meet */
373: cells[0] = cl; cells[1] = cr;
374: DMPlexGetMeet(*dm, 2, cells, &meetSize, &meet);
375: if (meetSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %d", f);
376: DMSetLabelValue(*dm, "Face Sets", meet[0], faceZoneIDs[f]);
377: DMPlexRestoreMeet(*dm, numFaceVertices, fverts, &meetSize, &meet);
378: } else {
379: for (fi = 0; fi < numFaceVertices; fi++) fverts[fi] = faces[f*numFaceEntries + fi] + numCells - 1;
380: DMPlexGetFullJoin(*dm, numFaceVertices, fverts, &joinSize, &join);
381: if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %d", f);
382: DMSetLabelValue(*dm, "Face Sets", join[0], faceZoneIDs[f]);
383: DMPlexRestoreJoin(*dm, numFaceVertices, fverts, &joinSize, &join);
384: }
385: }
386: PetscFree(fverts);
387: }
389: /* Read coordinates */
390: DMGetCoordinateSection(*dm, &coordSection);
391: PetscSectionSetNumFields(coordSection, 1);
392: PetscSectionSetFieldComponents(coordSection, 0, dim);
393: PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
394: for (v = numCells; v < numCells+numVertices; ++v) {
395: PetscSectionSetDof(coordSection, v, dim);
396: PetscSectionSetFieldDof(coordSection, v, 0, dim);
397: }
398: PetscSectionSetUp(coordSection);
399: PetscSectionGetStorageSize(coordSection, &coordSize);
400: VecCreate(PETSC_COMM_SELF, &coordinates);
401: PetscObjectSetName((PetscObject) coordinates, "coordinates");
402: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
403: VecSetType(coordinates, VECSTANDARD);
404: VecGetArray(coordinates, &coords);
405: if (!rank && coordsIn) {
406: for (v = 0; v < numVertices; ++v) {
407: for (d = 0; d < dim; ++d) {
408: coords[v*dim+d] = coordsIn[v*dim+d];
409: }
410: }
411: }
412: VecRestoreArray(coordinates, &coords);
413: DMSetCoordinatesLocal(*dm, coordinates);
414: VecDestroy(&coordinates);
415: if (!rank) {
416: if (cellVertices) {PetscFree(cellVertices);}
417: PetscFree(faces);
418: PetscFree(faceZoneIDs);
419: PetscFree(coordsIn);
420: }
421: return(0);
422: }