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: }