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

 10:   Input Parameters:
 11: + comm        - The MPI communicator
 12: . filename    - Name of the Fluent mesh file
 13: - interpolate - Create faces and edges in the mesh

 15:   Output Parameter:
 16: . dm  - The `DM` object representing the mesh

 18:   Level: beginner

 20: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexCreateFromFile()`, `DMPlexCreateFluent()`, `DMPlexCreate()`
 21: @*/
 22: PetscErrorCode DMPlexCreateFluentFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
 23: {
 24:   PetscViewer viewer;

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

 36: static PetscErrorCode DMPlexCreateFluent_ReadString(PetscViewer viewer, char *buffer, char delim)
 37: {
 38:   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:   if (!ret) buffer[i - 1] = '\0';
 43:   else 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;

 53:   if (binary) {
 54:     /* Extract raw file descriptor to read binary block */
 55:     PetscViewerASCIIGetPointer(viewer, &file);
 56:     fflush(file);
 57:     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);
 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, NULL, 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, NULL, 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;

 99:   /* Fast-forward to next section and derive its index */
100:   DMPlexCreateFluent_ReadString(viewer, buffer, '(');
101:   DMPlexCreateFluent_ReadString(viewer, buffer, ' ');
102:   snum = sscanf(buffer, "%d", &(s->index));
103:   /* If we can't match an index return -1 to signal end-of-file */
104:   if (snum < 1) {
105:     s->index = -1;
106:     return 0;
107:   }

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));

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));
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));
134:     if (s->zoneID == 0) { /* Header section */
135:       snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
137:     } else { /* Data section */
138:       snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
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));
156:     if (s->zoneID == 0) { /* Header section */
157:       snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
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));
163:       DMPlexCreateFluent_ReadString(viewer, buffer, '(');
164:       switch (s->nd) {
165:       case 0:
166:         numEntries = PETSC_DETERMINE;
167:         break;
168:       case 2:
169:         numEntries = 2 + 2;
170:         break; /* linear */
171:       case 3:
172:         numEntries = 2 + 3;
173:         break; /* triangular */
174:       case 4:
175:         numEntries = 2 + 4;
176:         break; /* quadrilateral */
177:       default:
178:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file");
179:       }
180:       numFaces = s->last - s->first + 1;
181:       if (numEntries != PETSC_DETERMINE) {
182:         /* Allocate space only if we already know the size of the block */
183:         PetscMalloc1(numEntries * numFaces, (PetscInt **)&s->data);
184:       }
185:       for (f = 0; f < numFaces; f++) {
186:         if (s->nd == 0) {
187:           /* Determine the size of the block for "mixed" facets */
188:           PetscInt numFaceVert = 0;
189:           DMPlexCreateFluent_ReadValues(viewer, &numFaceVert, 1, PETSC_INT, s->index == 2013 ? PETSC_TRUE : PETSC_FALSE);
190:           if (numEntries == PETSC_DETERMINE) {
191:             numEntries = numFaceVert + 2;
192:             PetscMalloc1(numEntries * numFaces, (PetscInt **)&s->data);
193:           } else {
195:           }
196:         }
197:         DMPlexCreateFluent_ReadValues(viewer, &(((PetscInt *)s->data)[f * numEntries]), numEntries, PETSC_INT, s->index == 2013 ? PETSC_TRUE : PETSC_FALSE);
198:       }
199:       s->nd = numEntries - 2;
200:       DMPlexCreateFluent_ReadString(viewer, buffer, ')');
201:     }
202:     DMPlexCreateFluent_ReadString(viewer, buffer, ')');

204:   } else { /* Unknown section type */
205:     PetscInt depth = 1;
206:     do {
207:       /* Match parentheses when parsing unknown sections */
208:       do PetscViewerRead(viewer, &(buffer[0]), 1, NULL, PETSC_CHAR);
209:       while (buffer[0] != '(' && buffer[0] != ')');
210:       if (buffer[0] == '(') depth++;
211:       if (buffer[0] == ')') depth--;
212:     } while (depth > 0);
213:     DMPlexCreateFluent_ReadString(viewer, buffer, '\n');
214:   }
215:   return 0;
216: }

218: /*@C
219:   DMPlexCreateFluent - Create a `DMPLEX` mesh from a Fluent mesh file.

221:   Collective

223:   Input Parameters:
224: + comm  - The MPI communicator
225: . viewer - The `PetscViewer` associated with a Fluent mesh file
226: - interpolate - Create faces and edges in the mesh

228:   Output Parameter:
229: . dm  - The `DM` object representing the mesh

231:   Note:
232:   http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm

234:   Level: beginner

236: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMCreate()`
237: @*/
238: PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
239: {
240:   PetscMPIInt  rank;
241:   PetscInt     c, v, dim = PETSC_DETERMINE, numCells = 0, numVertices = 0, numCellVertices = PETSC_DETERMINE;
242:   PetscInt     numFaces = PETSC_DETERMINE, f, numFaceEntries = PETSC_DETERMINE, numFaceVertices = PETSC_DETERMINE;
243:   PetscInt    *faces = NULL, *cellVertices = NULL, *faceZoneIDs = NULL;
244:   DMLabel      faceSets = NULL;
245:   PetscInt     d, coordSize;
246:   PetscScalar *coords, *coordsIn = NULL;
247:   PetscSection coordSection;
248:   Vec          coordinates;

250:   MPI_Comm_rank(comm, &rank);

252:   if (rank == 0) {
253:     FluentSection s;
254:     do {
255:       DMPlexCreateFluent_ReadSection(viewer, &s);
256:       if (s.index == 2) { /* Dimension */
257:         dim = s.nd;

259:       } else if (s.index == 10 || s.index == 2010) { /* Vertices */
260:         if (s.zoneID == 0) numVertices = s.last;
261:         else {
263:           coordsIn = (PetscScalar *)s.data;
264:         }

266:       } else if (s.index == 12 || s.index == 2012) { /* Cells */
267:         if (s.zoneID == 0) numCells = s.last;
268:         else {
269:           switch (s.nd) {
270:           case 0:
271:             numCellVertices = PETSC_DETERMINE;
272:             break;
273:           case 1:
274:             numCellVertices = 3;
275:             break; /* triangular */
276:           case 2:
277:             numCellVertices = 4;
278:             break; /* tetrahedral */
279:           case 3:
280:             numCellVertices = 4;
281:             break; /* quadrilateral */
282:           case 4:
283:             numCellVertices = 8;
284:             break; /* hexahedral */
285:           case 5:
286:             numCellVertices = 5;
287:             break; /* pyramid */
288:           case 6:
289:             numCellVertices = 6;
290:             break; /* wedge */
291:           default:
292:             numCellVertices = PETSC_DETERMINE;
293:           }
294:         }

296:       } else if (s.index == 13 || s.index == 2013) { /* Facets */
297:         if (s.zoneID == 0) {                         /* Header section */
298:           numFaces = (PetscInt)(s.last - s.first + 1);
299:           if (s.nd == 0 || s.nd == 5) numFaceVertices = PETSC_DETERMINE;
300:           else numFaceVertices = s.nd;
301:         } else { /* Data section */
302:           unsigned int z;

306:           if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd;
307:           numFaceEntries = numFaceVertices + 2;
308:           if (!faces) PetscMalloc1(numFaces * numFaceEntries, &faces);
309:           if (!faceZoneIDs) PetscMalloc1(numFaces, &faceZoneIDs);
310:           PetscMemcpy(&faces[(s.first - 1) * numFaceEntries], s.data, (s.last - s.first + 1) * numFaceEntries * sizeof(PetscInt));
311:           /* Record the zoneID for each face set */
312:           for (z = s.first - 1; z < s.last; z++) faceZoneIDs[z] = s.zoneID;
313:           PetscFree(s.data);
314:         }
315:       }
316:     } while (s.index >= 0);
317:   }
318:   MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);

321:   /* Allocate cell-vertex mesh */
322:   DMCreate(comm, dm);
323:   DMSetType(*dm, DMPLEX);
324:   DMSetDimension(*dm, dim);
325:   DMPlexSetChart(*dm, 0, numCells + numVertices);
326:   if (rank == 0) {
328:     /* If no cell type was given we assume simplices */
329:     if (numCellVertices == PETSC_DETERMINE) numCellVertices = numFaceVertices + 1;
330:     for (c = 0; c < numCells; ++c) DMPlexSetConeSize(*dm, c, numCellVertices);
331:   }
332:   DMSetUp(*dm);

334:   if (rank == 0 && faces) {
335:     /* Derive cell-vertex list from face-vertex and face-cell maps */
336:     PetscMalloc1(numCells * numCellVertices, &cellVertices);
337:     for (c = 0; c < numCells * numCellVertices; c++) cellVertices[c] = -1;
338:     for (f = 0; f < numFaces; f++) {
339:       PetscInt       *cell;
340:       const PetscInt  cl   = faces[f * numFaceEntries + numFaceVertices];
341:       const PetscInt  cr   = faces[f * numFaceEntries + numFaceVertices + 1];
342:       const PetscInt *face = &(faces[f * numFaceEntries]);

344:       if (cl > 0) {
345:         cell = &(cellVertices[(cl - 1) * numCellVertices]);
346:         for (v = 0; v < numFaceVertices; v++) {
347:           PetscBool found = PETSC_FALSE;
348:           for (c = 0; c < numCellVertices; c++) {
349:             if (cell[c] < 0) break;
350:             if (cell[c] == face[v] - 1 + numCells) {
351:               found = PETSC_TRUE;
352:               break;
353:             }
354:           }
355:           if (!found) cell[c] = face[v] - 1 + numCells;
356:         }
357:       }
358:       if (cr > 0) {
359:         cell = &(cellVertices[(cr - 1) * numCellVertices]);
360:         for (v = 0; v < numFaceVertices; v++) {
361:           PetscBool found = PETSC_FALSE;
362:           for (c = 0; c < numCellVertices; c++) {
363:             if (cell[c] < 0) break;
364:             if (cell[c] == face[v] - 1 + numCells) {
365:               found = PETSC_TRUE;
366:               break;
367:             }
368:           }
369:           if (!found) cell[c] = face[v] - 1 + numCells;
370:         }
371:       }
372:     }
373:     for (c = 0; c < numCells; c++) DMPlexSetCone(*dm, c, &(cellVertices[c * numCellVertices]));
374:   }
375:   DMPlexSymmetrize(*dm);
376:   DMPlexStratify(*dm);
377:   if (interpolate) {
378:     DM idm;

380:     DMPlexInterpolate(*dm, &idm);
381:     DMDestroy(dm);
382:     *dm = idm;
383:   }

385:   if (rank == 0 && faces) {
386:     PetscInt        fi, joinSize, meetSize, *fverts, cells[2];
387:     const PetscInt *join, *meet;
388:     PetscMalloc1(numFaceVertices, &fverts);
389:     /* Mark facets by finding the full join of all adjacent vertices */
390:     for (f = 0; f < numFaces; f++) {
391:       const PetscInt cl = faces[f * numFaceEntries + numFaceVertices] - 1;
392:       const PetscInt cr = faces[f * numFaceEntries + numFaceVertices + 1] - 1;
393:       if (cl > 0 && cr > 0) {
394:         /* If we know both adjoining cells we can use a single-level meet */
395:         cells[0] = cl;
396:         cells[1] = cr;
397:         DMPlexGetMeet(*dm, 2, cells, &meetSize, &meet);
399:         DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", meet[0], faceZoneIDs[f]);
400:         DMPlexRestoreMeet(*dm, numFaceVertices, fverts, &meetSize, &meet);
401:       } else {
402:         for (fi = 0; fi < numFaceVertices; fi++) fverts[fi] = faces[f * numFaceEntries + fi] + numCells - 1;
403:         DMPlexGetFullJoin(*dm, numFaceVertices, fverts, &joinSize, &join);
405:         DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], faceZoneIDs[f]);
406:         DMPlexRestoreJoin(*dm, numFaceVertices, fverts, &joinSize, &join);
407:       }
408:     }
409:     PetscFree(fverts);
410:   }

412:   { /* Create Face Sets label at all processes */
413:     enum {
414:       n = 1
415:     };
416:     PetscBool flag[n];

418:     flag[0] = faceSets ? PETSC_TRUE : PETSC_FALSE;
419:     MPI_Bcast(flag, n, MPIU_BOOL, 0, comm);
420:     if (flag[0]) DMCreateLabel(*dm, "Face Sets");
421:   }

423:   /* Read coordinates */
424:   DMGetCoordinateSection(*dm, &coordSection);
425:   PetscSectionSetNumFields(coordSection, 1);
426:   PetscSectionSetFieldComponents(coordSection, 0, dim);
427:   PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
428:   for (v = numCells; v < numCells + numVertices; ++v) {
429:     PetscSectionSetDof(coordSection, v, dim);
430:     PetscSectionSetFieldDof(coordSection, v, 0, dim);
431:   }
432:   PetscSectionSetUp(coordSection);
433:   PetscSectionGetStorageSize(coordSection, &coordSize);
434:   VecCreate(PETSC_COMM_SELF, &coordinates);
435:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
436:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
437:   VecSetType(coordinates, VECSTANDARD);
438:   VecGetArray(coordinates, &coords);
439:   if (rank == 0 && coordsIn) {
440:     for (v = 0; v < numVertices; ++v) {
441:       for (d = 0; d < dim; ++d) coords[v * dim + d] = coordsIn[v * dim + d];
442:     }
443:   }
444:   VecRestoreArray(coordinates, &coords);
445:   DMSetCoordinatesLocal(*dm, coordinates);
446:   VecDestroy(&coordinates);

448:   if (rank == 0) {
449:     PetscFree(cellVertices);
450:     PetscFree(faces);
451:     PetscFree(faceZoneIDs);
452:     PetscFree(coordsIn);
453:   }
454:   return 0;
455: }