Actual source code: plexfluent.c

petsc-3.9.4 2018-09-11
Report Typos and Errors
  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: }