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'; else buffer[i] = '\0';
 40:   return 0;
 41: }

 43: static PetscErrorCode DMPlexCreateFluent_ReadValues(PetscViewer viewer, void *data, PetscInt count, PetscDataType dtype, PetscBool binary)
 44: {
 45:   int            fdes=0;
 46:   FILE          *file;
 47:   PetscInt       i;

 49:   if (binary) {
 50:     /* Extract raw file descriptor to read binary block */
 51:     PetscViewerASCIIGetPointer(viewer, &file);
 52:     fflush(file); fdes = fileno(file);
 53:   }

 55:   if (!binary && dtype == PETSC_INT) {
 56:     char         cbuf[256];
 57:     unsigned int ibuf;
 58:     int          snum;
 59:     /* Parse hexadecimal ascii integers */
 60:     for (i = 0; i < count; i++) {
 61:       PetscViewerRead(viewer, cbuf, 1, NULL, PETSC_STRING);
 62:       snum = sscanf(cbuf, "%x", &ibuf);
 64:       ((PetscInt*)data)[i] = (PetscInt)ibuf;
 65:     }
 66:   } else if (binary && dtype == PETSC_INT) {
 67:     /* Always read 32-bit ints and cast to PetscInt */
 68:     int *ibuf;
 69:     PetscMalloc1(count, &ibuf);
 70:     PetscBinaryRead(fdes, ibuf, count, NULL, PETSC_ENUM);
 71:     PetscByteSwap(ibuf, PETSC_ENUM, count);
 72:     for (i = 0; i < count; i++) ((PetscInt*)data)[i] = (PetscInt)(ibuf[i]);
 73:     PetscFree(ibuf);

 75:  } else if (binary && dtype == PETSC_SCALAR) {
 76:     float *fbuf;
 77:     /* Always read 32-bit floats and cast to PetscScalar */
 78:     PetscMalloc1(count, &fbuf);
 79:     PetscBinaryRead(fdes, fbuf, count, NULL, PETSC_FLOAT);
 80:     PetscByteSwap(fbuf, PETSC_FLOAT, count);
 81:     for (i = 0; i < count; i++) ((PetscScalar*)data)[i] = (PetscScalar)(fbuf[i]);
 82:     PetscFree(fbuf);
 83:   } else {
 84:     PetscViewerASCIIRead(viewer, data, count, NULL, dtype);
 85:   }
 86:   return 0;
 87: }

 89: static PetscErrorCode DMPlexCreateFluent_ReadSection(PetscViewer viewer, FluentSection *s)
 90: {
 91:   char           buffer[PETSC_MAX_PATH_LEN];
 92:   int            snum;

 94:   /* Fast-forward to next section and derive its index */
 95:   DMPlexCreateFluent_ReadString(viewer, buffer, '(');
 96:   DMPlexCreateFluent_ReadString(viewer, buffer, ' ');
 97:   snum = sscanf(buffer, "%d", &(s->index));
 98:   /* If we can't match an index return -1 to signal end-of-file */
 99:   if (snum < 1) {s->index = -1;   return 0;}

101:   if (s->index == 0) {           /* Comment */
102:     DMPlexCreateFluent_ReadString(viewer, buffer, ')');

104:   } else if (s->index == 2) {    /* Dimension */
105:     DMPlexCreateFluent_ReadString(viewer, buffer, ')');
106:     snum = sscanf(buffer, "%d", &(s->nd));

109:   } else if (s->index == 10 || s->index == 2010) {   /* Vertices */
110:     DMPlexCreateFluent_ReadString(viewer, buffer, ')');
111:     snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
113:     if (s->zoneID > 0) {
114:       PetscInt numCoords = s->last - s->first + 1;
115:       DMPlexCreateFluent_ReadString(viewer, buffer, '(');
116:       PetscMalloc1(s->nd*numCoords, (PetscScalar**)&s->data);
117:       DMPlexCreateFluent_ReadValues(viewer, s->data, s->nd*numCoords, PETSC_SCALAR, s->index==2010 ? PETSC_TRUE : PETSC_FALSE);
118:       DMPlexCreateFluent_ReadString(viewer, buffer, ')');
119:     }
120:     DMPlexCreateFluent_ReadString(viewer, buffer, ')');

122:   } else if (s->index == 12 || s->index == 2012) {   /* Cells */
123:     DMPlexCreateFluent_ReadString(viewer, buffer, ')');
124:     snum = sscanf(buffer, "(%x", &(s->zoneID));
126:     if (s->zoneID == 0) {  /* Header section */
127:       snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
129:     } else {               /* Data section */
130:       snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
132:       if (s->nd == 0) {
133:         /* Read cell type definitions for mixed cells */
134:         PetscInt numCells = s->last - s->first + 1;
135:         DMPlexCreateFluent_ReadString(viewer, buffer, '(');
136:         PetscMalloc1(numCells, (PetscInt**)&s->data);
137:         DMPlexCreateFluent_ReadValues(viewer, s->data, numCells, PETSC_INT, s->index==2012 ? PETSC_TRUE : PETSC_FALSE);
138:         PetscFree(s->data);
139:         DMPlexCreateFluent_ReadString(viewer, buffer, ')');
140:       }
141:     }
142:     DMPlexCreateFluent_ReadString(viewer, buffer, ')');

144:   } else if (s->index == 13 || s->index == 2013) {   /* Faces */
145:     DMPlexCreateFluent_ReadString(viewer, buffer, ')');
146:     snum = sscanf(buffer, "(%x", &(s->zoneID));
148:     if (s->zoneID == 0) {  /* Header section */
149:       snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
151:     } else {               /* Data section */
152:       PetscInt f, numEntries, numFaces;
153:       snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
155:       DMPlexCreateFluent_ReadString(viewer, buffer, '(');
156:       switch (s->nd) {
157:       case 0: numEntries = PETSC_DETERMINE; break;
158:       case 2: numEntries = 2 + 2; break;  /* linear */
159:       case 3: numEntries = 2 + 3; break;  /* triangular */
160:       case 4: numEntries = 2 + 4; break;  /* quadrilateral */
161:       default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file");
162:       }
163:       numFaces = s->last-s->first + 1;
164:       if (numEntries != PETSC_DETERMINE) {
165:         /* Allocate space only if we already know the size of the block */
166:         PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);
167:       }
168:       for (f = 0; f < numFaces; f++) {
169:         if (s->nd == 0) {
170:           /* Determine the size of the block for "mixed" facets */
171:           PetscInt numFaceVert = 0;
172:           DMPlexCreateFluent_ReadValues(viewer, &numFaceVert, 1, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);
173:           if (numEntries == PETSC_DETERMINE) {
174:             numEntries = numFaceVert + 2;
175:             PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);
176:           } else {
178:           }
179:         }
180:         DMPlexCreateFluent_ReadValues(viewer, &(((PetscInt*)s->data)[f*numEntries]), numEntries, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);
181:       }
182:       s->nd = numEntries - 2;
183:       DMPlexCreateFluent_ReadString(viewer, buffer, ')');
184:     }
185:     DMPlexCreateFluent_ReadString(viewer, buffer, ')');

187:   } else {                       /* Unknown section type */
188:     PetscInt depth = 1;
189:     do {
190:       /* Match parentheses when parsing unknown sections */
191:       do PetscViewerRead(viewer, &(buffer[0]), 1, NULL, PETSC_CHAR);
192:       while (buffer[0] != '(' && buffer[0] != ')');
193:       if (buffer[0] == '(') depth++;
194:       if (buffer[0] == ')') depth--;
195:     } while (depth > 0);
196:     DMPlexCreateFluent_ReadString(viewer, buffer, '\n');
197:   }
198:   return 0;
199: }

201: /*@C
202:   DMPlexCreateFluent - Create a DMPlex mesh from a Fluent mesh file.

204:   Collective

206:   Input Parameters:
207: + comm  - The MPI communicator
208: . viewer - The Viewer associated with a Fluent mesh file
209: - interpolate - Create faces and edges in the mesh

211:   Output Parameter:
212: . dm  - The DM object representing the mesh

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

216:   Level: beginner

218: .seealso: DMPLEX, DMCreate()
219: @*/
220: PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
221: {
222:   PetscMPIInt    rank;
223:   PetscInt       c, v, dim = PETSC_DETERMINE, numCells = 0, numVertices = 0, numCellVertices = PETSC_DETERMINE;
224:   PetscInt       numFaces = PETSC_DETERMINE, f, numFaceEntries = PETSC_DETERMINE, numFaceVertices = PETSC_DETERMINE;
225:   PetscInt      *faces = NULL, *cellVertices = NULL, *faceZoneIDs = NULL;
226:   DMLabel        faceSets = NULL;
227:   PetscInt       d, coordSize;
228:   PetscScalar   *coords, *coordsIn = NULL;
229:   PetscSection   coordSection;
230:   Vec            coordinates;

232:   MPI_Comm_rank(comm, &rank);

234:   if (rank == 0) {
235:     FluentSection s;
236:     do {
237:       DMPlexCreateFluent_ReadSection(viewer, &s);
238:       if (s.index == 2) {                 /* Dimension */
239:         dim = s.nd;

241:       } else if (s.index == 10 || s.index == 2010) { /* Vertices */
242:         if (s.zoneID == 0) numVertices = s.last;
243:         else {
245:           coordsIn = (PetscScalar *) s.data;
246:         }

248:       } else if (s.index == 12 || s.index == 2012) { /* Cells */
249:         if (s.zoneID == 0) numCells = s.last;
250:         else {
251:           switch (s.nd) {
252:           case 0: numCellVertices = PETSC_DETERMINE; break;
253:           case 1: numCellVertices = 3; break;  /* triangular */
254:           case 2: numCellVertices = 4; break;  /* tetrahedral */
255:           case 3: numCellVertices = 4; break;  /* quadrilateral */
256:           case 4: numCellVertices = 8; break;  /* hexahedral */
257:           case 5: numCellVertices = 5; break;  /* pyramid */
258:           case 6: numCellVertices = 6; break;  /* wedge */
259:           default: numCellVertices = PETSC_DETERMINE;
260:           }
261:         }

263:       } else if (s.index == 13 || s.index == 2013) { /* Facets */
264:         if (s.zoneID == 0) {  /* Header section */
265:           numFaces = (PetscInt) (s.last - s.first + 1);
266:           if (s.nd == 0 || s.nd == 5) numFaceVertices = PETSC_DETERMINE;
267:           else numFaceVertices = s.nd;
268:         } else {              /* Data section */
269:           unsigned int z;

273:           if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd;
274:           numFaceEntries = numFaceVertices + 2;
275:           if (!faces) PetscMalloc1(numFaces*numFaceEntries, &faces);
276:           if (!faceZoneIDs) PetscMalloc1(numFaces, &faceZoneIDs);
277:           PetscMemcpy(&faces[(s.first-1)*numFaceEntries], s.data, (s.last-s.first+1)*numFaceEntries*sizeof(PetscInt));
278:           /* Record the zoneID for each face set */
279:           for (z = s.first -1; z < s.last; z++) faceZoneIDs[z] = s.zoneID;
280:           PetscFree(s.data);
281:         }
282:       }
283:     } while (s.index >= 0);
284:   }
285:   MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);

288:   /* Allocate cell-vertex mesh */
289:   DMCreate(comm, dm);
290:   DMSetType(*dm, DMPLEX);
291:   DMSetDimension(*dm, dim);
292:   DMPlexSetChart(*dm, 0, numCells + numVertices);
293:   if (rank == 0) {
295:     /* If no cell type was given we assume simplices */
296:     if (numCellVertices == PETSC_DETERMINE) numCellVertices = numFaceVertices + 1;
297:     for (c = 0; c < numCells; ++c) DMPlexSetConeSize(*dm, c, numCellVertices);
298:   }
299:   DMSetUp(*dm);

301:   if (rank == 0 && faces) {
302:     /* Derive cell-vertex list from face-vertex and face-cell maps */
303:     PetscMalloc1(numCells*numCellVertices, &cellVertices);
304:     for (c = 0; c < numCells*numCellVertices; c++) cellVertices[c] = -1;
305:     for (f = 0; f < numFaces; f++) {
306:       PetscInt *cell;
307:       const PetscInt cl = faces[f*numFaceEntries + numFaceVertices];
308:       const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1];
309:       const PetscInt *face = &(faces[f*numFaceEntries]);

311:       if (cl > 0) {
312:         cell = &(cellVertices[(cl-1) * numCellVertices]);
313:         for (v = 0; v < numFaceVertices; v++) {
314:           PetscBool found = PETSC_FALSE;
315:           for (c = 0; c < numCellVertices; c++) {
316:             if (cell[c] < 0) break;
317:             if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
318:           }
319:           if (!found) cell[c] = face[v]-1 + numCells;
320:         }
321:       }
322:       if (cr > 0) {
323:         cell = &(cellVertices[(cr-1) * numCellVertices]);
324:         for (v = 0; v < numFaceVertices; v++) {
325:           PetscBool found = PETSC_FALSE;
326:           for (c = 0; c < numCellVertices; c++) {
327:             if (cell[c] < 0) break;
328:             if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
329:           }
330:           if (!found) cell[c] = face[v]-1 + numCells;
331:         }
332:       }
333:     }
334:     for (c = 0; c < numCells; c++) {
335:       DMPlexSetCone(*dm, c, &(cellVertices[c*numCellVertices]));
336:     }
337:   }
338:   DMPlexSymmetrize(*dm);
339:   DMPlexStratify(*dm);
340:   if (interpolate) {
341:     DM idm;

343:     DMPlexInterpolate(*dm, &idm);
344:     DMDestroy(dm);
345:     *dm  = idm;
346:   }

348:   if (rank == 0 && faces) {
349:     PetscInt fi, joinSize, meetSize, *fverts, cells[2];
350:     const PetscInt *join, *meet;
351:     PetscMalloc1(numFaceVertices, &fverts);
352:     /* Mark facets by finding the full join of all adjacent vertices */
353:     for (f = 0; f < numFaces; f++) {
354:       const PetscInt cl = faces[f*numFaceEntries + numFaceVertices] - 1;
355:       const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1] - 1;
356:       if (cl > 0 && cr > 0) {
357:         /* If we know both adjoining cells we can use a single-level meet */
358:         cells[0] = cl; cells[1] = cr;
359:         DMPlexGetMeet(*dm, 2, cells, &meetSize, &meet);
361:         DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", meet[0], faceZoneIDs[f]);
362:         DMPlexRestoreMeet(*dm, numFaceVertices, fverts, &meetSize, &meet);
363:       } else {
364:         for (fi = 0; fi < numFaceVertices; fi++) fverts[fi] = faces[f*numFaceEntries + fi] + numCells - 1;
365:         DMPlexGetFullJoin(*dm, numFaceVertices, fverts, &joinSize, &join);
367:         DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], faceZoneIDs[f]);
368:         DMPlexRestoreJoin(*dm, numFaceVertices, fverts, &joinSize, &join);
369:       }
370:     }
371:     PetscFree(fverts);
372:   }

374:   { /* Create Face Sets label at all processes */
375:     enum {n = 1};
376:     PetscBool flag[n];

378:     flag[0] = faceSets ? PETSC_TRUE : PETSC_FALSE;
379:     MPI_Bcast(flag, n, MPIU_BOOL, 0, comm);
380:     if (flag[0]) DMCreateLabel(*dm, "Face Sets");
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 == 0 && 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);

410:   if (rank == 0) {
411:     PetscFree(cellVertices);
412:     PetscFree(faces);
413:     PetscFree(faceZoneIDs);
414:     PetscFree(coordsIn);
415:   }
416:   return 0;
417: }