Actual source code: plexfluent.c

  1: #define PETSC_DESIRE_FEATURE_TEST_MACROS /* for fileno() */
  2: #include <petsc/private/dmpleximpl.h>

  4: /* Utility struct to store the contents of a Fluent file in memory */
  5: typedef struct {
  6:   int          index; /* Type of section */
  7:   unsigned int zoneID;
  8:   unsigned int first;
  9:   unsigned int last;
 10:   int          type;
 11:   int          nd; /* Either ND or element-type */
 12:   void        *data;
 13: } FluentSection;

 15: /*@
 16:   DMPlexCreateFluentFromFile - Create a `DMPLEX` mesh from a Fluent mesh file

 18:   Collective

 20:   Input Parameters:
 21: + comm        - The MPI communicator
 22: . filename    - Name of the Fluent mesh file
 23: - interpolate - Create faces and edges in the mesh

 25:   Output Parameter:
 26: . dm - The `DM` object representing the mesh

 28:   Level: beginner

 30: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateFromFile()`, `DMPlexCreateFluent()`, `DMPlexCreate()`
 31: @*/
 32: PetscErrorCode DMPlexCreateFluentFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
 33: {
 34:   PetscViewer viewer;

 36:   PetscFunctionBegin;
 37:   /* Create file viewer and build plex */
 38:   PetscCall(PetscViewerCreate(comm, &viewer));
 39:   PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII));
 40:   PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ));
 41:   PetscCall(PetscViewerFileSetName(viewer, filename));
 42:   PetscCall(DMPlexCreateFluent(comm, viewer, interpolate, dm));
 43:   PetscCall(PetscViewerDestroy(&viewer));
 44:   PetscFunctionReturn(PETSC_SUCCESS);
 45: }

 47: static PetscErrorCode DMPlexCreateFluent_ReadString(PetscViewer viewer, char *buffer, char delim)
 48: {
 49:   PetscInt ret, i = 0;

 51:   PetscFunctionBegin;
 52:   do PetscCall(PetscViewerRead(viewer, &buffer[i++], 1, &ret, PETSC_CHAR));
 53:   while (ret > 0 && buffer[i - 1] != '\0' && buffer[i - 1] != delim && i < PETSC_MAX_PATH_LEN - 1);
 54:   if (!ret) buffer[i - 1] = '\0';
 55:   else buffer[i] = '\0';
 56:   PetscCheck(i < PETSC_MAX_PATH_LEN - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Buffer overflow! This is not a valid Fluent file.");
 57:   PetscFunctionReturn(PETSC_SUCCESS);
 58: }

 60: static PetscErrorCode DMPlexCreateFluent_ReadValues(PetscViewer viewer, void *data, PetscInt count, PetscDataType dtype, PetscBool binary, PetscInt *numClosingParens)
 61: {
 62:   int      fdes = 0;
 63:   FILE    *file;
 64:   PetscInt i;

 66:   PetscFunctionBegin;
 67:   *numClosingParens = 0;
 68:   if (binary) {
 69:     /* Extract raw file descriptor to read binary block */
 70:     PetscCall(PetscViewerASCIIGetPointer(viewer, &file));
 71:     PetscCall(PetscFFlush(file));
 72:     fdes = fileno(file);
 73:   }

 75:   if (!binary && dtype == PETSC_INT) {
 76:     char         cbuf[256];
 77:     unsigned int ibuf;
 78:     int          snum;
 79:     /* Parse hexadecimal ascii integers */
 80:     for (i = 0; i < count; i++) {
 81:       size_t len;

 83:       PetscCall(PetscViewerRead(viewer, cbuf, 1, NULL, PETSC_STRING));
 84:       snum = sscanf(cbuf, "%x", &ibuf);
 85:       PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
 86:       ((PetscInt *)data)[i] = (PetscInt)ibuf;
 87:       // Check for trailing parentheses
 88:       PetscCall(PetscStrlen(cbuf, &len));
 89:       while (cbuf[len - 1] == ')' && len > 0) {
 90:         ++(*numClosingParens);
 91:         --len;
 92:       }
 93:     }
 94:   } else if (binary && dtype == PETSC_INT) {
 95:     /* Always read 32-bit ints and cast to PetscInt */
 96:     int *ibuf;
 97:     PetscCall(PetscMalloc1(count, &ibuf));
 98:     PetscCall(PetscBinaryRead(fdes, ibuf, count, NULL, PETSC_ENUM));
 99:     PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, count));
100:     for (i = 0; i < count; i++) ((PetscInt *)data)[i] = ibuf[i];
101:     PetscCall(PetscFree(ibuf));

103:   } else if (binary && dtype == PETSC_SCALAR) {
104:     float *fbuf;
105:     /* Always read 32-bit floats and cast to PetscScalar */
106:     PetscCall(PetscMalloc1(count, &fbuf));
107:     PetscCall(PetscBinaryRead(fdes, fbuf, count, NULL, PETSC_FLOAT));
108:     PetscCall(PetscByteSwap(fbuf, PETSC_FLOAT, count));
109:     for (i = 0; i < count; i++) ((PetscScalar *)data)[i] = fbuf[i];
110:     PetscCall(PetscFree(fbuf));
111:   } else {
112:     PetscCall(PetscViewerASCIIRead(viewer, data, count, NULL, dtype));
113:   }
114:   PetscFunctionReturn(PETSC_SUCCESS);
115: }

117: static PetscErrorCode DMPlexCreateFluent_ReadSection(PetscViewer viewer, FluentSection *s)
118: {
119:   char buffer[PETSC_MAX_PATH_LEN];
120:   int  snum;

122:   PetscFunctionBegin;
123:   /* Fast-forward to next section and derive its index */
124:   PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
125:   PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ' '));
126:   snum = sscanf(buffer, "%d", &s->index);
127:   /* If we can't match an index return -1 to signal end-of-file */
128:   if (snum < 1) {
129:     s->index = -1;
130:     PetscFunctionReturn(PETSC_SUCCESS);
131:   }

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

136:   } else if (s->index == 2) { /* Dimension */
137:     PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
138:     snum = sscanf(buffer, "%d", &s->nd);
139:     PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");

141:   } else if (s->index == 10 || s->index == 2010) { /* Vertices */
142:     PetscInt numClosingParens = 0;

144:     PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
145:     snum = sscanf(buffer, "(%x %x %x %d %d)", &s->zoneID, &s->first, &s->last, &s->type, &s->nd);
146:     PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
147:     if (s->zoneID > 0) {
148:       PetscInt numCoords = s->last - s->first + 1;
149:       PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
150:       PetscCall(PetscMalloc1(s->nd * numCoords, (PetscScalar **)&s->data));
151:       PetscCall(DMPlexCreateFluent_ReadValues(viewer, s->data, s->nd * numCoords, PETSC_SCALAR, s->index == 2010 ? PETSC_TRUE : PETSC_FALSE, &numClosingParens));
152:       if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
153:       else --numClosingParens;
154:     }
155:     if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
156:     else --numClosingParens;
157:     PetscCheck(!numClosingParens, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
158:   } else if (s->index == 12 || s->index == 2012) { /* Cells */
159:     PetscInt numClosingParens = 0;

161:     PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
162:     snum = sscanf(buffer, "(%x", &s->zoneID);
163:     PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
164:     if (s->zoneID == 0) { /* Header section */
165:       snum = sscanf(buffer, "(%x %x %x %d)", &s->zoneID, &s->first, &s->last, &s->nd);
166:       PetscCheck(snum == 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
167:     } else { /* Data section */
168:       snum = sscanf(buffer, "(%x %x %x %d %d)", &s->zoneID, &s->first, &s->last, &s->type, &s->nd);
169:       PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
170:       if (s->nd == 0) {
171:         /* Read cell type definitions for mixed cells */
172:         PetscInt numCells = s->last - s->first + 1;
173:         PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
174:         PetscCall(PetscMalloc1(numCells, (PetscInt **)&s->data));
175:         PetscCall(DMPlexCreateFluent_ReadValues(viewer, s->data, numCells, PETSC_INT, s->index == 2012 ? PETSC_TRUE : PETSC_FALSE, &numClosingParens));
176:         if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
177:         else --numClosingParens;
178:       }
179:     }
180:     if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
181:     else --numClosingParens;
182:     PetscCheck(!numClosingParens, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
183:   } else if (s->index == 13 || s->index == 2013) { /* Faces */
184:     PetscInt numClosingParens = 0;

186:     PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
187:     snum = sscanf(buffer, "(%x", &s->zoneID);
188:     PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
189:     if (s->zoneID == 0) { /* Header section */
190:       snum = sscanf(buffer, "(%x %x %x %d)", &s->zoneID, &s->first, &s->last, &s->nd);
191:       PetscCheck(snum == 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
192:     } else { /* Data section */
193:       PetscInt numEntries, numFaces, maxsize = 0, offset = 0;

195:       snum = sscanf(buffer, "(%x %x %x %d %d)", &s->zoneID, &s->first, &s->last, &s->type, &s->nd);
196:       PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
197:       PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
198:       switch (s->nd) {
199:       case 0:
200:         numEntries = PETSC_DETERMINE;
201:         break;
202:       case 2:
203:         numEntries = 2 + 2;
204:         break; /* linear */
205:       case 3:
206:         numEntries = 2 + 3;
207:         break; /* triangular */
208:       case 4:
209:         numEntries = 2 + 4;
210:         break; /* quadrilateral */
211:       default:
212:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file");
213:       }
214:       numFaces = s->last - s->first + 1;
215:       if (numEntries != PETSC_DETERMINE) {
216:         /* Allocate space only if we already know the size of the block */
217:         PetscCall(PetscMalloc1(numEntries * numFaces, (PetscInt **)&s->data));
218:       }
219:       for (PetscInt f = 0; f < numFaces; f++) {
220:         if (s->nd == 0) {
221:           /* Determine the size of the block for "mixed" facets */
222:           PetscInt numFaceVert = 0;
223:           PetscCall(DMPlexCreateFluent_ReadValues(viewer, &numFaceVert, 1, PETSC_INT, s->index == 2013 ? PETSC_TRUE : PETSC_FALSE, &numClosingParens));
224:           if (!f) {
225:             maxsize = (numFaceVert + 3) * numFaces;
226:             PetscCall(PetscMalloc1(maxsize, (PetscInt **)&s->data));
227:           } else {
228:             if (offset + numFaceVert + 3 >= maxsize) {
229:               PetscInt *tmp;

231:               PetscCall(PetscMalloc1(maxsize * 2, &tmp));
232:               PetscCall(PetscArraycpy(tmp, (PetscInt *)s->data, maxsize));
233:               PetscCall(PetscFree(s->data));
234:               maxsize *= 2;
235:               s->data = tmp;
236:             }
237:           }
238:           ((PetscInt *)s->data)[offset] = numFaceVert;
239:           ++offset;
240:           numEntries = numFaceVert + 2;
241:         }
242:         PetscCall(DMPlexCreateFluent_ReadValues(viewer, &(((PetscInt *)s->data)[offset]), numEntries, PETSC_INT, s->index == 2013 ? PETSC_TRUE : PETSC_FALSE, &numClosingParens));
243:         offset += numEntries;
244:       }
245:       if (s->nd != 0) PetscCall(PetscMPIIntCast(numEntries - 2, &s->nd));
246:       if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
247:       else --numClosingParens;
248:     }
249:     if (!numClosingParens) PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
250:     else --numClosingParens;
251:     PetscCheck(!numClosingParens, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
252:   } else if (s->index == 39) { /* Label information */
253:     char labelName[PETSC_MAX_PATH_LEN];
254:     char caseName[PETSC_MAX_PATH_LEN];

256:     PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
257:     snum = sscanf(buffer, "(%u %s %s %d)", &s->zoneID, caseName, labelName, &s->nd);
258:     PetscCheck(snum == 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file: %d", snum);
259:     PetscInt depth = 1;
260:     do {
261:       /* Match parentheses when parsing unknown sections */
262:       do PetscCall(PetscViewerRead(viewer, &buffer[0], 1, NULL, PETSC_CHAR));
263:       while (buffer[0] != '(' && buffer[0] != ')');
264:       if (buffer[0] == '(') depth++;
265:       if (buffer[0] == ')') depth--;
266:     } while (depth > 0);
267:     PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '\n'));
268:     PetscCall(PetscStrallocpy(labelName, (char **)&s->data));
269:     PetscCall(PetscInfo((PetscObject)viewer, "CASE: Zone ID %u is label %s\n", s->zoneID, labelName));
270:   } else { /* Unknown section type */
271:     PetscInt depth = 1;
272:     do {
273:       /* Match parentheses when parsing unknown sections */
274:       do PetscCall(PetscViewerRead(viewer, &buffer[0], 1, NULL, PETSC_CHAR));
275:       while (buffer[0] != '(' && buffer[0] != ')');
276:       if (buffer[0] == '(') depth++;
277:       if (buffer[0] == ')') depth--;
278:     } while (depth > 0);
279:     PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '\n'));
280:   }
281:   PetscFunctionReturn(PETSC_SUCCESS);
282: }

284: // Inserts point `face` with orientation `ornt` into the cone of point `cell` at position `c`, which is the first empty slot
285: static PetscErrorCode InsertFace(DM dm, PetscInt cell, PetscInt face, PetscInt ornt)
286: {
287:   const PetscInt *cone;
288:   PetscInt        coneSize, c;

290:   PetscFunctionBegin;
291:   PetscCall(DMPlexGetCone(dm, cell, &cone));
292:   PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
293:   for (c = 0; c < coneSize; ++c)
294:     if (cone[c] < 0) break;
295:   PetscCheck(c < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " could not be inserted in cone of cell %" PetscInt_FMT " with size %" PetscInt_FMT, face, cell, coneSize);
296:   PetscCall(DMPlexInsertCone(dm, cell, c, face));
297:   PetscCall(DMPlexInsertConeOrientation(dm, cell, c, ornt));
298:   PetscFunctionReturn(PETSC_SUCCESS);
299: }

301: static PetscErrorCode ReorderPolygon(DM dm, PetscInt cell)
302: {
303:   const PetscInt *cone, *ornt;
304:   PetscInt        coneSize, newCone[16], newOrnt[16];

306:   PetscFunctionBegin;
307:   PetscCall(DMPlexGetOrientedCone(dm, cell, &cone, &ornt));
308:   PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
309:   newCone[0] = cone[0];
310:   newOrnt[0] = ornt[0];
311:   for (PetscInt c = 1; c < coneSize; ++c) {
312:     const PetscInt *fcone;
313:     PetscInt        firstVertex, lastVertex, c2;

315:     PetscCall(DMPlexGetCone(dm, newCone[c - 1], &fcone));
316:     lastVertex = newOrnt[c - 1] ? fcone[0] : fcone[1];
317:     for (c2 = 0; c2 < coneSize; ++c2) {
318:       const PetscInt *fcone2;

320:       PetscCall(DMPlexGetCone(dm, cone[c2], &fcone2));
321:       firstVertex = ornt[c2] ? fcone2[1] : fcone2[0];
322:       if (lastVertex == firstVertex) {
323:         // Point `cell` matched point `lastVertex` on face `cone[c2]` with orientation `ornt[c2]`
324:         break;
325:       }
326:     }
327:     PetscCheck(c2 < coneSize, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " could not find a face match as position %" PetscInt_FMT, cell, c);
328:     newCone[c] = cone[c2];
329:     newOrnt[c] = ornt[c2];
330:   }
331:   {
332:     const PetscInt *fcone, *fcone2;
333:     PetscInt        vertex, vertex2;

335:     PetscCall(DMPlexGetCone(dm, newCone[coneSize - 1], &fcone));
336:     PetscCall(DMPlexGetCone(dm, newCone[0], &fcone2));
337:     vertex  = newOrnt[coneSize - 1] ? fcone[0] : fcone[1];
338:     vertex2 = newOrnt[0] ? fcone2[1] : fcone2[0];
339:     PetscCheck(vertex == vertex2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " did not match at the endpoint", cell);
340:   }
341:   PetscCall(DMPlexSetCone(dm, cell, newCone));
342:   PetscCall(DMPlexSetConeOrientation(dm, cell, newOrnt));
343:   PetscCall(DMPlexRestoreOrientedCone(dm, cell, &cone, &ornt));
344:   PetscFunctionReturn(PETSC_SUCCESS);
345: }

347: static PetscErrorCode ReorderTetrahedron(PetscViewer viewer, DM dm, PetscInt cell)
348: {
349:   const PetscInt *cone, *ornt, *fcone, *fornt, *farr, faces[4] = {0, 1, 3, 2};
350:   PetscInt        newCone[16], newOrnt[16];

352:   PetscFunctionBegin;
353:   PetscCall(DMPlexGetOrientedCone(dm, cell, &cone, &ornt));
354:   newCone[0] = cone[0];
355:   newOrnt[0] = ornt[0];
356:   PetscCall(DMPlexGetOrientedCone(dm, newCone[0], &fcone, &fornt));
357:   farr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, newOrnt[0]);
358:   // Loop over each edge in the initial triangle
359:   for (PetscInt e = 0; e < 3; ++e) {
360:     const PetscInt edge = fcone[farr[e * 2 + 0]], eornt = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[e * 2 + 1], fornt[farr[e * 2 + 0]]);
361:     PetscInt       c;

363:     // Loop over each remaining face in the tetrahedron
364:     //   On face `newCone[0]`, trying to match edge `edge` with final orientation `eornt` to an edge on another face
365:     for (c = 1; c < 4; ++c) {
366:       const PetscInt *fcone2, *fornt2, *farr2;
367:       PetscInt        c2;
368:       PetscBool       flip = PETSC_FALSE;

370:       // Checking face `cone[c]` with orientation `ornt[c]`
371:       PetscCall(DMPlexGetOrientedCone(dm, cone[c], &fcone2, &fornt2));
372:       farr2 = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, ornt[c]);
373:       // Check for edge
374:       for (c2 = 0; c2 < 3; ++c2) {
375:         const PetscInt edge2 = fcone2[farr2[c2 * 2 + 0]], eornt2 = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr2[c2 * 2 + 1], fornt2[farr2[c2 * 2 + 0]]);
376:         // Trying to match edge `edge2` with final orientation `eornt2`
377:         if (edge == edge2) {
378:           //PetscCheck(eornt == -(eornt2 + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell % " PetscInt_FMT " edge %" PetscInt_FMT " (%" PetscInt_FMT ") found twice with the same orientation in face %" PetscInt_FMT " edge %" PetscInt_FMT, cell, edge, e, c, c2);
379:           // Matched face `newCone[0]` with orientation `newOrnt[0]` to face `cone[c]` with orientation `ornt[c]` along edge `edge`
380:           PetscCall(PetscInfo((PetscObject)viewer, "CASE: Matched cell %" PetscInt_FMT " edge %" PetscInt_FMT "/%" PetscInt_FMT " (%" PetscInt_FMT ") to face %" PetscInt_FMT "/%" PetscInt_FMT " edge %" PetscInt_FMT " (%" PetscInt_FMT ")\n", cell, edge, e, eornt, cone[c], c, c2, eornt2));
381:           flip = eornt != -(eornt2 + 1) ? PETSC_TRUE : PETSC_FALSE;
382:           break;
383:         }
384:       }
385:       if (c2 < 3) {
386:         newCone[faces[e + 1]] = cone[c];
387:         // Compute new orientation of face based on which edge was matched (only the first edge matches a side different from 0)
388:         //   Face 1 should match its edge 2
389:         //   Face 2 should match its edge 0
390:         //   Face 3 should match its edge 0
391:         if (flip) {
392:           newOrnt[faces[e + 1]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_TRIANGLE, -((c2 + (!e ? 1 : 2)) % 3 + 1), ornt[c]);
393:         } else {
394:           newOrnt[faces[e + 1]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_TRIANGLE, !e ? (c2 + 1) % 3 : c2, ornt[c]);
395:         }
396:         break;
397:       }
398:     }
399:     PetscCheck(c < 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " could not find a face match for edge %" PetscInt_FMT, cell, e);
400:   }
401:   PetscCall(DMPlexRestoreOrientedCone(dm, newCone[0], &fcone, &fornt));
402:   PetscCall(DMPlexSetCone(dm, cell, newCone));
403:   PetscCall(DMPlexSetConeOrientation(dm, cell, newOrnt));
404:   PetscCall(DMPlexRestoreOrientedCone(dm, cell, &cone, &ornt));
405:   PetscFunctionReturn(PETSC_SUCCESS);
406: }

408: static PetscErrorCode ReorderHexahedron(DM dm, PetscInt cell)
409: {
410:   const PetscInt *cone, *ornt, *fcone, *fornt, *farr;
411:   const PetscInt  faces[6] = {0, 5, 3, 4, 2, 1};
412:   PetscInt        used[6]  = {1, 0, 0, 0, 0, 0};
413:   PetscInt        newCone[16], newOrnt[16];

415:   PetscFunctionBegin;
416:   PetscCall(DMPlexGetOrientedCone(dm, cell, &cone, &ornt));
417:   newCone[0] = cone[0];
418:   newOrnt[0] = ornt[0];
419:   PetscCall(DMPlexGetOrientedCone(dm, newCone[0], &fcone, &fornt));
420:   farr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, newOrnt[0]);
421:   // Loop over each edge in the initial quadrilateral
422:   for (PetscInt e = 0; e < 4; ++e) {
423:     const PetscInt edge = fcone[farr[e * 2 + 0]], eornt = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[e * 2 + 1], fornt[farr[e * 2 + 0]]);
424:     PetscInt       c;

426:     // Loop over each remaining face in the hexahedron
427:     //   On face `newCone[0]`, trying to match edge `edge` with final orientation `eornt` to an edge on another face
428:     for (c = 1; c < 6; ++c) {
429:       const PetscInt *fcone2, *fornt2, *farr2;
430:       PetscInt        c2;

432:       // Checking face `cone[c]` with orientation `ornt[c]`
433:       PetscCall(DMPlexGetOrientedCone(dm, cone[c], &fcone2, &fornt2));
434:       farr2 = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, ornt[c]);
435:       // Check for edge
436:       for (c2 = 0; c2 < 4; ++c2) {
437:         const PetscInt edge2 = fcone2[farr2[c2 * 2 + 0]], eornt2 = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr2[c2 * 2 + 1], fornt2[farr2[c2 * 2 + 0]]);
438:         // Trying to match edge `edge2` with final orientation `eornt2`
439:         if (edge == edge2) {
440:           PetscCheck(eornt == -(eornt2 + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " found twice with the same orientation", edge);
441:           // Matched face `newCone[0]` with orientation `newOrnt[0]` to face `cone[c]` with orientation `ornt[c]` along edge `edge`
442:           break;
443:         }
444:       }
445:       if (c2 < 4) {
446:         used[c]               = 1;
447:         newCone[faces[e + 1]] = cone[c];
448:         // Compute new orientation of face based on which edge was matched (only the first edge matches a side different from 0)
449:         newOrnt[faces[e + 1]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_QUADRILATERAL, !e ? (c2 + 1) % 4 : c2, ornt[c]);
450:         break;
451:       }
452:     }
453:     PetscCheck(c < 6, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " could not find a face match for edge %" PetscInt_FMT, cell, e);
454:   }
455:   PetscCall(DMPlexRestoreOrientedCone(dm, newCone[0], &fcone, &fornt));
456:   // Add last face
457:   {
458:     PetscInt c, c2;

460:     for (c = 1; c < 6; ++c)
461:       if (!used[c]) break;
462:     PetscCheck(c < 6, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " could not find an available face", cell);
463:     // Match first edge to 3rd edge in newCone[2]
464:     {
465:       const PetscInt *fcone2, *fornt2, *farr2;

467:       PetscCall(DMPlexGetOrientedCone(dm, newCone[2], &fcone, &fornt));
468:       farr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, newOrnt[2]);
469:       PetscCall(DMPlexGetOrientedCone(dm, cone[c], &fcone2, &fornt2));
470:       farr2 = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, ornt[c]);

472:       const PetscInt e    = 2;
473:       const PetscInt edge = fcone[farr[e * 2 + 0]], eornt = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[e * 2 + 1], fornt[farr[e * 2 + 0]]);
474:       // Trying to match edge `edge` with final orientation `eornt` of face `newCone[2]` to some edge of face `cone[c]` with orientation `ornt[c]`
475:       for (c2 = 0; c2 < 4; ++c2) {
476:         const PetscInt edge2 = fcone2[farr2[c2 * 2 + 0]], eornt2 = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr2[c2 * 2 + 1], fornt2[farr2[c2 * 2 + 0]]);
477:         // Trying to match edge `edge2` with final orientation `eornt2`
478:         if (edge == edge2) {
479:           PetscCheck(eornt == -(eornt2 + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " found twice with the same orientation", edge);
480:           // Matched face `newCone[2]` with orientation `newOrnt[2]` to face `cone[c]` with orientation `ornt[c]` along edge `edge`
481:           break;
482:         }
483:       }
484:       PetscCheck(c2 < 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not fit last face in");
485:     }
486:     newCone[faces[5]] = cone[c];
487:     // Compute new orientation of face based on which edge was matched
488:     newOrnt[faces[5]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_QUADRILATERAL, c2, ornt[c]);
489:     PetscCall(DMPlexRestoreOrientedCone(dm, newCone[0], &fcone, &fornt));
490:   }
491:   PetscCall(DMPlexSetCone(dm, cell, newCone));
492:   PetscCall(DMPlexSetConeOrientation(dm, cell, newOrnt));
493:   PetscCall(DMPlexRestoreOrientedCone(dm, cell, &cone, &ornt));
494:   PetscFunctionReturn(PETSC_SUCCESS);
495: }

497: // {0, 1, 2}, {3, 4, 5}, {0, 2, 4, 3}, {2, 1, 5, 4}, {1, 0, 3, 5}
498: static PetscErrorCode ReorderWedge(DM dm, PetscInt cell)
499: {
500:   const PetscInt *cone, *ornt, *fcone, *fornt, *farr;
501:   const PetscInt  faces[5] = {0, 4, 3, 2, 1};
502:   PetscInt        used[5]  = {0, 0, 0, 0, 0};
503:   PetscInt        newCone[16], newOrnt[16], cS, bottom = 0;

505:   PetscFunctionBegin;
506:   PetscCall(DMPlexGetConeSize(dm, cell, &cS));
507:   PetscCall(DMPlexGetOrientedCone(dm, cell, &cone, &ornt));
508:   for (PetscInt c = 0; c < cS; ++c) {
509:     DMPolytopeType ct;

511:     PetscCall(DMPlexGetCellType(dm, cone[c], &ct));
512:     if (ct == DM_POLYTOPE_TRIANGLE) {
513:       bottom = c;
514:       break;
515:     }
516:   }
517:   used[bottom] = 1;
518:   newCone[0]   = cone[bottom];
519:   newOrnt[0]   = ornt[bottom];
520:   PetscCall(DMPlexGetOrientedCone(dm, newCone[0], &fcone, &fornt));
521:   farr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, newOrnt[0]);
522:   // Loop over each edge in the initial triangle
523:   for (PetscInt e = 0; e < 3; ++e) {
524:     const PetscInt edge = fcone[farr[e * 2 + 0]], eornt = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[e * 2 + 1], fornt[farr[e * 2 + 0]]);
525:     PetscInt       c;

527:     // Loop over each remaining face in the prism
528:     //   On face `newCone[0]`, trying to match edge `edge` with final orientation `eornt` to an edge on another face
529:     for (c = 0; c < 5; ++c) {
530:       const PetscInt *fcone2, *fornt2, *farr2;
531:       DMPolytopeType  ct;
532:       PetscInt        c2;

534:       if (c == bottom) continue;
535:       PetscCall(DMPlexGetCellType(dm, cone[c], &ct));
536:       if (ct != DM_POLYTOPE_QUADRILATERAL) continue;
537:       // Checking face `cone[c]` with orientation `ornt[c]`
538:       PetscCall(DMPlexGetOrientedCone(dm, cone[c], &fcone2, &fornt2));
539:       farr2 = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, ornt[c]);
540:       // Check for edge
541:       for (c2 = 0; c2 < 4; ++c2) {
542:         const PetscInt edge2 = fcone2[farr2[c2 * 2 + 0]], eornt2 = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr2[c2 * 2 + 1], fornt2[farr2[c2 * 2 + 0]]);
543:         // Trying to match edge `edge2` with final orientation `eornt2`
544:         if (edge == edge2) {
545:           PetscCheck(eornt == -(eornt2 + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " found twice with the same orientation", edge);
546:           // Matched face `newCone[0]` with orientation `newOrnt[0]` to face `cone[c]` with orientation `ornt[c]` along edge `edge`
547:           break;
548:         }
549:       }
550:       if (c2 < 4) {
551:         used[c]               = 1;
552:         newCone[faces[e + 1]] = cone[c];
553:         // Compute new orientation of face based on which edge was matched, edge 0 should always match the bottom
554:         newOrnt[faces[e + 1]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_QUADRILATERAL, c2, ornt[c]);
555:         break;
556:       }
557:     }
558:     PetscCheck(c < 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " could not find a face match for edge %" PetscInt_FMT, cell, e);
559:   }
560:   PetscCall(DMPlexRestoreOrientedCone(dm, newCone[0], &fcone, &fornt));
561:   // Add last face
562:   {
563:     PetscInt c, c2;

565:     for (c = 0; c < 5; ++c)
566:       if (!used[c]) break;
567:     PetscCheck(c < 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell %" PetscInt_FMT " could not find an available face", cell);
568:     // Match first edge to 3rd edge in newCone[2]
569:     {
570:       const PetscInt *fcone2, *fornt2, *farr2;

572:       PetscCall(DMPlexGetOrientedCone(dm, newCone[2], &fcone, &fornt));
573:       farr = DMPolytopeTypeGetArrangement(DM_POLYTOPE_QUADRILATERAL, newOrnt[2]);
574:       PetscCall(DMPlexGetOrientedCone(dm, cone[c], &fcone2, &fornt2));
575:       farr2 = DMPolytopeTypeGetArrangement(DM_POLYTOPE_TRIANGLE, ornt[c]);

577:       const PetscInt e    = 2;
578:       const PetscInt edge = fcone[farr[e * 2 + 0]], eornt = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr[e * 2 + 1], fornt[farr[e * 2 + 0]]);
579:       // Trying to match edge `edge` with final orientation `eornt` of face `newCone[2]` to some edge of face `cone[c]` with orientation `ornt[c]`
580:       for (c2 = 0; c2 < 3; ++c2) {
581:         const PetscInt edge2 = fcone2[farr2[c2 * 2 + 0]], eornt2 = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_SEGMENT, farr2[c2 * 2 + 1], fornt2[farr2[c2 * 2 + 0]]);
582:         // Trying to match edge `edge2` with final orientation `eornt2`
583:         if (edge == edge2) {
584:           PetscCheck(eornt == -(eornt2 + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %" PetscInt_FMT " found twice with the same orientation", edge);
585:           // Matched face `newCone[2]` with orientation `newOrnt[2]` to face `cone[c]` with orientation `ornt[c]` along edge `edge`
586:           break;
587:         }
588:       }
589:       PetscCheck(c2 < 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not fit last face in");
590:     }
591:     newCone[faces[4]] = cone[c];
592:     // Compute new orientation of face based on which edge was matched
593:     newOrnt[faces[4]] = DMPolytopeTypeComposeOrientation(DM_POLYTOPE_TRIANGLE, c2, ornt[c]);
594:     PetscCall(DMPlexRestoreOrientedCone(dm, newCone[0], &fcone, &fornt));
595:   }
596:   PetscCall(DMPlexSetCone(dm, cell, newCone));
597:   PetscCall(DMPlexSetConeOrientation(dm, cell, newOrnt));
598:   PetscCall(DMPlexRestoreOrientedCone(dm, cell, &cone, &ornt));
599:   PetscFunctionReturn(PETSC_SUCCESS);
600: }

602: static PetscErrorCode ReorderCell(PetscViewer viewer, DM dm, PetscInt cell, DMPolytopeType ct)
603: {
604:   PetscFunctionBegin;
605:   switch (ct) {
606:   case DM_POLYTOPE_TRIANGLE:
607:   case DM_POLYTOPE_QUADRILATERAL:
608:     PetscCall(ReorderPolygon(dm, cell));
609:     break;
610:   case DM_POLYTOPE_TETRAHEDRON:
611:     PetscCall(ReorderTetrahedron(viewer, dm, cell));
612:     break;
613:   case DM_POLYTOPE_HEXAHEDRON:
614:     PetscCall(ReorderHexahedron(dm, cell));
615:     break;
616:   case DM_POLYTOPE_TRI_PRISM:
617:     PetscCall(ReorderWedge(dm, cell));
618:     break;
619:   default:
620:     PetscCheck(0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Celltype %s is unsupported", DMPolytopeTypes[ct]);
621:     break;
622:   }
623:   PetscFunctionReturn(PETSC_SUCCESS);
624: }

626: static PetscErrorCode GetNumCellFaces(int nd, PetscInt *numCellFaces, DMPolytopeType *ct)
627: {
628:   PetscFunctionBegin;
629:   *ct = DM_POLYTOPE_POINT;
630:   switch (nd) {
631:   case 0:
632:     *numCellFaces = PETSC_DETERMINE;
633:     break;
634:   case 1:
635:     *numCellFaces = 3;
636:     *ct           = DM_POLYTOPE_TRIANGLE;
637:     break;
638:   case 2:
639:     *numCellFaces = 4;
640:     *ct           = DM_POLYTOPE_TETRAHEDRON;
641:     break;
642:   case 3:
643:     *numCellFaces = 4;
644:     *ct           = DM_POLYTOPE_QUADRILATERAL;
645:     break;
646:   case 4:
647:     *numCellFaces = 6;
648:     *ct           = DM_POLYTOPE_HEXAHEDRON;
649:     break;
650:   case 5:
651:     *numCellFaces = 5;
652:     *ct           = DM_POLYTOPE_PYRAMID;
653:     break;
654:   case 6:
655:     *numCellFaces = 5;
656:     *ct           = DM_POLYTOPE_TRI_PRISM;
657:     break;
658:   default:
659:     *numCellFaces = PETSC_DETERMINE;
660:   }
661:   PetscFunctionReturn(PETSC_SUCCESS);
662: }

664: /*@C
665:   DMPlexCreateFluent - Create a `DMPLEX` mesh from a Fluent mesh file <http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm>.

667:   Collective

669:   Input Parameters:
670: + comm        - The MPI communicator
671: . viewer      - The `PetscViewer` associated with a Fluent mesh file
672: - interpolate - Create faces and edges in the mesh

674:   Output Parameter:
675: . dm - The `DM` object representing the mesh

677:   Level: beginner

679: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`
680: @*/
681: PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
682: {
683:   PetscInt        dim          = PETSC_DETERMINE;
684:   PetscInt        numCells     = 0;
685:   PetscInt        numVertices  = 0;
686:   PetscInt       *cellSizes    = NULL;
687:   DMPolytopeType *cellTypes    = NULL;
688:   PetscInt        numFaces     = 0;
689:   PetscInt       *faces        = NULL;
690:   PetscInt       *faceSizes    = NULL;
691:   PetscInt       *faceAdjCell  = NULL;
692:   PetscInt       *cellVertices = NULL;
693:   unsigned int   *faceZoneIDs  = NULL;
694:   DMLabel         faceSets     = NULL;
695:   DMLabel        *zoneLabels   = NULL;
696:   const char    **zoneNames    = NULL;
697:   unsigned int    maxZoneID    = 0;
698:   PetscScalar    *coordsIn     = NULL;
699:   PetscScalar    *coords;
700:   PetscSection    coordSection;
701:   Vec             coordinates;
702:   PetscInt        coordSize, maxFaceSize = 0, totFaceVert = 0, f;
703:   PetscMPIInt     rank;

705:   PetscFunctionBegin;
706:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

708:   if (rank == 0) {
709:     FluentSection s;

711:     s.data   = NULL;
712:     numFaces = PETSC_DETERMINE;
713:     do {
714:       PetscCall(DMPlexCreateFluent_ReadSection(viewer, &s));
715:       if (s.index == 2) { /* Dimension */
716:         dim = s.nd;
717:         PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found dimension: %" PetscInt_FMT "\n", dim));
718:       } else if (s.index == 10 || s.index == 2010) { /* Vertices */
719:         if (s.zoneID == 0) {
720:           numVertices = s.last;
721:           PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found number of vertices: %" PetscInt_FMT "\n", numVertices));
722:         } else {
723:           PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found vertex coordinates\n"));
724:           PetscCheck(!coordsIn, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files");
725:           coordsIn = (PetscScalar *)s.data;
726:         }

728:       } else if (s.index == 12 || s.index == 2012) { /* Cells */
729:         if (s.zoneID == 0) {
730:           numCells = s.last;
731:           PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found number of cells %" PetscInt_FMT "\n", numCells));
732:         } else {
733:           PetscCall(PetscMalloc2(numCells, &cellSizes, numCells, &cellTypes));
734:           for (PetscInt c = 0; c < numCells; ++c) PetscCall(GetNumCellFaces(s.nd ? s.nd : (int)((PetscInt *)s.data)[c], &cellSizes[c], &cellTypes[c]));
735:           PetscCall(PetscFree(s.data));
736:           PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found number of cell faces %" PetscInt_FMT "\n", numCells && s.nd ? cellSizes[0] : 0));
737:         }
738:       } else if (s.index == 13 || s.index == 2013) { /* Facets */
739:         if (s.zoneID == 0) {                         /* Header section */
740:           numFaces = (PetscInt)(s.last - s.first + 1);
741:           PetscCall(PetscInfo((PetscObject)viewer, "CASE: Found number of faces %" PetscInt_FMT " face vertices: %d\n", numFaces, s.nd));
742:         } else { /* Data section */
743:           PetscInt *tmp;
744:           PetscInt  totSize = 0, offset = 0, doffset;

746:           PetscCheck(numFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file");
747:           if (!faceZoneIDs) PetscCall(PetscMalloc3(numFaces, &faceSizes, numFaces * 2, &faceAdjCell, numFaces, &faceZoneIDs));
748:           // Record the zoneID and face size for each face set
749:           for (unsigned int z = s.first - 1; z < s.last; z++) {
750:             faceZoneIDs[z] = s.zoneID;
751:             if (s.nd) {
752:               faceSizes[z] = s.nd;
753:             } else {
754:               faceSizes[z] = ((PetscInt *)s.data)[offset];
755:               offset += faceSizes[z] + 3;
756:             }
757:             totSize += faceSizes[z];
758:             maxFaceSize = PetscMax(maxFaceSize, faceSizes[z]);
759:           }

761:           offset  = totFaceVert;
762:           doffset = s.nd ? 0 : 1;
763:           PetscCall(PetscMalloc1(totFaceVert + totSize, &tmp));
764:           if (faces) PetscCall(PetscArraycpy(tmp, faces, totFaceVert));
765:           PetscCall(PetscFree(faces));
766:           totFaceVert += totSize;
767:           faces = tmp;

769:           // Record face vertices and adjacent faces
770:           const PetscInt Nfz = s.last - s.first + 1;
771:           for (PetscInt f = 0; f < Nfz; ++f) {
772:             const PetscInt face     = f + s.first - 1;
773:             const PetscInt faceSize = faceSizes[face];

775:             for (PetscInt v = 0; v < faceSize; ++v) faces[offset + v] = ((PetscInt *)s.data)[doffset + v];
776:             faceAdjCell[face * 2 + 0] = ((PetscInt *)s.data)[doffset + faceSize + 0];
777:             faceAdjCell[face * 2 + 1] = ((PetscInt *)s.data)[doffset + faceSize + 1];
778:             offset += faceSize;
779:             doffset += faceSize + (s.nd ? 2 : 3);
780:           }
781:           PetscCall(PetscFree(s.data));
782:         }
783:       } else if (s.index == 39) { /* Label information */
784:         if (s.zoneID >= maxZoneID) {
785:           DMLabel     *tmpL;
786:           const char **tmp;
787:           unsigned int newmax = maxZoneID + 1;

789:           while (newmax < s.zoneID + 1) newmax *= 2;
790:           PetscCall(PetscCalloc2(newmax, &tmp, newmax, &tmpL));
791:           for (PetscInt i = 0; i < (PetscInt)maxZoneID; ++i) {
792:             tmp[i]  = zoneNames[i];
793:             tmpL[i] = zoneLabels[i];
794:           }
795:           maxZoneID = newmax;
796:           PetscCall(PetscFree2(zoneNames, zoneLabels));
797:           zoneNames  = tmp;
798:           zoneLabels = tmpL;
799:         }
800:         zoneNames[s.zoneID] = (const char *)s.data;
801:       }
802:     } while (s.index >= 0);
803:   }
804:   PetscCallMPI(MPI_Bcast(&dim, 1, MPIU_INT, 0, comm));
805:   PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension");

807:   /* Allocate cell-vertex mesh */
808:   PetscCall(DMCreate(comm, dm));
809:   PetscCall(DMSetType(*dm, DMPLEX));
810:   PetscCall(DMSetDimension(*dm, dim));
811:   // We do not want this label automatically computed, instead we fill it here
812:   PetscCall(DMCreateLabel(*dm, "celltype"));
813:   PetscCall(DMPlexSetChart(*dm, 0, numCells + numFaces + numVertices));
814:   if (rank == 0) {
815:     PetscCheck(numCells >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown number of cells in Fluent file");
816:     for (PetscInt c = 0; c < numCells; ++c) {
817:       PetscCall(DMPlexSetConeSize(*dm, c, cellSizes[c]));
818:       PetscCall(DMPlexSetCellType(*dm, c, cellTypes[c]));
819:     }
820:     for (PetscInt v = numCells; v < numCells + numVertices; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT));
821:     for (PetscInt f = 0; f < numFaces; ++f) {
822:       DMPolytopeType ct;

824:       switch (faceSizes[f]) {
825:       case 2:
826:         ct = DM_POLYTOPE_SEGMENT;
827:         break;
828:       case 3:
829:         ct = DM_POLYTOPE_TRIANGLE;
830:         break;
831:       case 4:
832:         ct = DM_POLYTOPE_QUADRILATERAL;
833:         break;
834:       default:
835:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file with cone size %" PetscInt_FMT, faceSizes[f]);
836:       }
837:       PetscCall(DMPlexSetConeSize(*dm, f + numCells + numVertices, faceSizes[f]));
838:       PetscCall(DMPlexSetCellType(*dm, f + numCells + numVertices, ct));
839:     }
840:   }
841:   PetscCall(DMSetUp(*dm));

843:   if (rank == 0 && faces) {
844:     PetscSection s;
845:     PetscInt    *cones, csize, foffset = 0;

847:     PetscCall(DMPlexGetCones(*dm, &cones));
848:     PetscCall(DMPlexGetConeSection(*dm, &s));
849:     PetscCall(PetscSectionGetConstrainedStorageSize(s, &csize));
850:     for (PetscInt c = 0; c < csize; ++c) cones[c] = -1;
851:     for (PetscInt f = 0; f < numFaces; f++) {
852:       const PetscInt cl   = faceAdjCell[f * 2 + 0] - 1;
853:       const PetscInt cr   = faceAdjCell[f * 2 + 1] - 1;
854:       const PetscInt face = f + numCells + numVertices;
855:       PetscInt       fcone[16];

857:       // How could Fluent define the outward normal differently? Is there no end to the pain?
858:       if (dim == 3) {
859:         if (cl >= 0) PetscCall(InsertFace(*dm, cl, face, -1));
860:         if (cr >= 0) PetscCall(InsertFace(*dm, cr, face, 0));
861:       } else {
862:         if (cl >= 0) PetscCall(InsertFace(*dm, cl, face, 0));
863:         if (cr >= 0) PetscCall(InsertFace(*dm, cr, face, -1));
864:       }
865:       PetscCheck(faceSizes[f] < 16, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of face vertices %" PetscInt_FMT " exceeds temporary storage", faceSizes[f]);
866:       for (PetscInt v = 0; v < faceSizes[f]; ++v) fcone[v] = faces[foffset + v] + numCells - 1;
867:       foffset += faceSizes[f];
868:       PetscCall(DMPlexSetCone(*dm, face, fcone));
869:     }
870:   }
871:   PetscCall(DMPlexSymmetrize(*dm));
872:   PetscCall(DMPlexStratify(*dm));
873:   if (dim == 3) {
874:     DM idm;

876:     PetscCall(DMCreate(PetscObjectComm((PetscObject)*dm), &idm));
877:     PetscCall(DMSetType(idm, DMPLEX));
878:     PetscCall(DMSetDimension(idm, dim));
879:     PetscCall(DMPlexInterpolateFaces_Internal(*dm, 1, idm));
880:     PetscCall(DMDestroy(dm));
881:     *dm = idm;
882:   }
883:   PetscCall(DMViewFromOptions(*dm, NULL, "-cas_dm_view"));
884:   if (rank == 0 && faces) {
885:     for (PetscInt c = 0; c < numCells; ++c) PetscCall(ReorderCell(viewer, *dm, c, cellTypes[c]));
886:   }

888:   if (rank == 0 && faces) {
889:     PetscInt        joinSize, meetSize, *fverts, cells[2];
890:     const PetscInt *join, *meet;
891:     PetscInt        foffset = 0;

893:     PetscCall(PetscMalloc1(maxFaceSize, &fverts));
894:     /* Mark facets by finding the full join of all adjacent vertices */
895:     for (f = 0; f < numFaces; f++) {
896:       const PetscInt cl = faceAdjCell[f * 2 + 0] - 1;
897:       const PetscInt cr = faceAdjCell[f * 2 + 1] - 1;
898:       const PetscInt id = (PetscInt)faceZoneIDs[f];

900:       if (cl > 0 && cr > 0) {
901:         /* If we know both adjoining cells we can use a single-level meet */
902:         cells[0] = cl;
903:         cells[1] = cr;
904:         PetscCall(DMPlexGetMeet(*dm, 2, cells, &meetSize, &meet));
905:         PetscCheck(meetSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %" PetscInt_FMT " cells: %" PetscInt_FMT ", %" PetscInt_FMT, f, cl, cr);
906:         PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", meet[0], id));
907:         if (zoneNames && zoneNames[id]) PetscCall(DMSetLabelValue_Fast(*dm, &zoneLabels[id], zoneNames[id], meet[0], 1));
908:         PetscCall(DMPlexRestoreMeet(*dm, meetSize, fverts, &meetSize, &meet));
909:       } else {
910:         for (PetscInt fi = 0; fi < faceSizes[f]; fi++) fverts[fi] = faces[foffset + fi] + numCells - 1;
911:         PetscCall(DMPlexGetFullJoin(*dm, faceSizes[f], fverts, &joinSize, &join));
912:         PetscCheck(joinSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %" PetscInt_FMT, f);
913:         PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], id));
914:         if (zoneNames && zoneNames[id]) PetscCall(DMSetLabelValue_Fast(*dm, &zoneLabels[id], zoneNames[id], join[0], 1));
915:         PetscCall(DMPlexRestoreJoin(*dm, joinSize, fverts, &joinSize, &join));
916:       }
917:       foffset += faceSizes[f];
918:     }
919:     PetscCall(PetscFree(fverts));
920:   }

922:   { /* Create Face Sets label at all processes */
923:     enum {
924:       n = 1
925:     };
926:     PetscBool flag[n];

928:     flag[0] = faceSets ? PETSC_TRUE : PETSC_FALSE;
929:     PetscCallMPI(MPI_Bcast(flag, n, MPI_C_BOOL, 0, comm));
930:     if (flag[0]) PetscCall(DMCreateLabel(*dm, "Face Sets"));
931:     // TODO Code to create all the zone labels on each process
932:   }

934:   if (!interpolate) {
935:     DM udm;

937:     PetscCall(DMPlexUninterpolate(*dm, &udm));
938:     PetscCall(DMDestroy(dm));
939:     *dm = udm;
940:   }

942:   /* Read coordinates */
943:   PetscCall(DMGetCoordinateSection(*dm, &coordSection));
944:   PetscCall(PetscSectionSetNumFields(coordSection, 1));
945:   PetscCall(PetscSectionSetFieldComponents(coordSection, 0, dim));
946:   PetscCall(PetscSectionSetChart(coordSection, numCells, numCells + numVertices));
947:   for (PetscInt v = numCells; v < numCells + numVertices; ++v) {
948:     PetscCall(PetscSectionSetDof(coordSection, v, dim));
949:     PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, dim));
950:   }
951:   PetscCall(PetscSectionSetUp(coordSection));
952:   PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize));
953:   PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates));
954:   PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
955:   PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE));
956:   PetscCall(VecSetType(coordinates, VECSTANDARD));
957:   PetscCall(VecGetArray(coordinates, &coords));
958:   if (rank == 0 && coordsIn) {
959:     for (PetscInt v = 0; v < numVertices; ++v) {
960:       for (PetscInt d = 0; d < dim; ++d) coords[v * dim + d] = coordsIn[v * dim + d];
961:     }
962:   }
963:   PetscCall(VecRestoreArray(coordinates, &coords));
964:   PetscCall(DMSetCoordinatesLocal(*dm, coordinates));
965:   PetscCall(VecDestroy(&coordinates));

967:   if (rank == 0) {
968:     PetscCall(PetscFree(cellVertices));
969:     PetscCall(PetscFree2(cellSizes, cellTypes));
970:     PetscCall(PetscFree(faces));
971:     PetscCall(PetscFree3(faceSizes, faceAdjCell, faceZoneIDs));
972:     PetscCall(PetscFree(coordsIn));
973:     if (zoneNames)
974:       for (PetscInt i = 0; i < (PetscInt)maxZoneID; ++i) PetscCall(PetscFree(zoneNames[i]));
975:     PetscCall(PetscFree2(zoneNames, zoneLabels));
976:   }
977:   PetscFunctionReturn(PETSC_SUCCESS);
978: }