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: /*@
      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: [](ch_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:   PetscFunctionBegin;
     27:   /* Create file viewer and build plex */
     28:   PetscCall(PetscViewerCreate(comm, &viewer));
     29:   PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII));
     30:   PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ));
     31:   PetscCall(PetscViewerFileSetName(viewer, filename));
     32:   PetscCall(DMPlexCreateFluent(comm, viewer, interpolate, dm));
     33:   PetscCall(PetscViewerDestroy(&viewer));
     34:   PetscFunctionReturn(PETSC_SUCCESS);
     35: }

     37: static PetscErrorCode DMPlexCreateFluent_ReadString(PetscViewer viewer, char *buffer, char delim)
     38: {
     39:   PetscInt ret, i = 0;

     41:   PetscFunctionBegin;
     42:   do PetscCall(PetscViewerRead(viewer, &buffer[i++], 1, &ret, PETSC_CHAR));
     43:   while (ret > 0 && buffer[i - 1] != '\0' && buffer[i - 1] != delim);
     44:   if (!ret) buffer[i - 1] = '\0';
     45:   else buffer[i] = '\0';
     46:   PetscFunctionReturn(PETSC_SUCCESS);
     47: }

     49: static PetscErrorCode DMPlexCreateFluent_ReadValues(PetscViewer viewer, void *data, PetscInt count, PetscDataType dtype, PetscBool binary)
     50: {
     51:   int      fdes = 0;
     52:   FILE    *file;
     53:   PetscInt i;

     55:   PetscFunctionBegin;
     56:   if (binary) {
     57:     /* Extract raw file descriptor to read binary block */
     58:     PetscCall(PetscViewerASCIIGetPointer(viewer, &file));
     59:     PetscCall(PetscFFlush(file));
     60:     fdes = fileno(file);
     61:   }

     63:   if (!binary && dtype == PETSC_INT) {
     64:     char         cbuf[256];
     65:     unsigned int ibuf;
     66:     int          snum;
     67:     /* Parse hexadecimal ascii integers */
     68:     for (i = 0; i < count; i++) {
     69:       PetscCall(PetscViewerRead(viewer, cbuf, 1, NULL, PETSC_STRING));
     70:       snum = sscanf(cbuf, "%x", &ibuf);
     71:       PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
     72:       ((PetscInt *)data)[i] = (PetscInt)ibuf;
     73:     }
     74:   } else if (binary && dtype == PETSC_INT) {
     75:     /* Always read 32-bit ints and cast to PetscInt */
     76:     int *ibuf;
     77:     PetscCall(PetscMalloc1(count, &ibuf));
     78:     PetscCall(PetscBinaryRead(fdes, ibuf, count, NULL, PETSC_ENUM));
     79:     PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, count));
     80:     for (i = 0; i < count; i++) ((PetscInt *)data)[i] = (PetscInt)ibuf[i];
     81:     PetscCall(PetscFree(ibuf));

     83:   } else if (binary && dtype == PETSC_SCALAR) {
     84:     float *fbuf;
     85:     /* Always read 32-bit floats and cast to PetscScalar */
     86:     PetscCall(PetscMalloc1(count, &fbuf));
     87:     PetscCall(PetscBinaryRead(fdes, fbuf, count, NULL, PETSC_FLOAT));
     88:     PetscCall(PetscByteSwap(fbuf, PETSC_FLOAT, count));
     89:     for (i = 0; i < count; i++) ((PetscScalar *)data)[i] = (PetscScalar)fbuf[i];
     90:     PetscCall(PetscFree(fbuf));
     91:   } else {
     92:     PetscCall(PetscViewerASCIIRead(viewer, data, count, NULL, dtype));
     93:   }
     94:   PetscFunctionReturn(PETSC_SUCCESS);
     95: }

     97: static PetscErrorCode DMPlexCreateFluent_ReadSection(PetscViewer viewer, FluentSection *s)
     98: {
     99:   char buffer[PETSC_MAX_PATH_LEN];
    100:   int  snum;

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

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

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

    121:   } else if (s->index == 10 || s->index == 2010) { /* Vertices */
    122:     PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
    123:     snum = sscanf(buffer, "(%x %x %x %d %d)", &s->zoneID, &s->first, &s->last, &s->type, &s->nd);
    124:     PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
    125:     if (s->zoneID > 0) {
    126:       PetscInt numCoords = s->last - s->first + 1;
    127:       PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
    128:       PetscCall(PetscMalloc1(s->nd * numCoords, (PetscScalar **)&s->data));
    129:       PetscCall(DMPlexCreateFluent_ReadValues(viewer, s->data, s->nd * numCoords, PETSC_SCALAR, s->index == 2010 ? PETSC_TRUE : PETSC_FALSE));
    130:       PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
    131:     }
    132:     PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));

    134:   } else if (s->index == 12 || s->index == 2012) { /* Cells */
    135:     PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
    136:     snum = sscanf(buffer, "(%x", &s->zoneID);
    137:     PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
    138:     if (s->zoneID == 0) { /* Header section */
    139:       snum = sscanf(buffer, "(%x %x %x %d)", &s->zoneID, &s->first, &s->last, &s->nd);
    140:       PetscCheck(snum == 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
    141:     } else { /* Data section */
    142:       snum = sscanf(buffer, "(%x %x %x %d %d)", &s->zoneID, &s->first, &s->last, &s->type, &s->nd);
    143:       PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
    144:       if (s->nd == 0) {
    145:         /* Read cell type definitions for mixed cells */
    146:         PetscInt numCells = s->last - s->first + 1;
    147:         PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
    148:         PetscCall(PetscMalloc1(numCells, (PetscInt **)&s->data));
    149:         PetscCall(DMPlexCreateFluent_ReadValues(viewer, s->data, numCells, PETSC_INT, s->index == 2012 ? PETSC_TRUE : PETSC_FALSE));
    150:         PetscCall(PetscFree(s->data));
    151:         PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
    152:       }
    153:     }
    154:     PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));

    156:   } else if (s->index == 13 || s->index == 2013) { /* Faces */
    157:     PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
    158:     snum = sscanf(buffer, "(%x", &s->zoneID);
    159:     PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
    160:     if (s->zoneID == 0) { /* Header section */
    161:       snum = sscanf(buffer, "(%x %x %x %d)", &s->zoneID, &s->first, &s->last, &s->nd);
    162:       PetscCheck(snum == 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
    163:     } else { /* Data section */
    164:       PetscInt f, numEntries, numFaces;
    165:       snum = sscanf(buffer, "(%x %x %x %d %d)", &s->zoneID, &s->first, &s->last, &s->type, &s->nd);
    166:       PetscCheck(snum == 5, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
    167:       PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, '('));
    168:       switch (s->nd) {
    169:       case 0:
    170:         numEntries = PETSC_DETERMINE;
    171:         break;
    172:       case 2:
    173:         numEntries = 2 + 2;
    174:         break; /* linear */
    175:       case 3:
    176:         numEntries = 2 + 3;
    177:         break; /* triangular */
    178:       case 4:
    179:         numEntries = 2 + 4;
    180:         break; /* quadrilateral */
    181:       default:
    182:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file");
    183:       }
    184:       numFaces = s->last - s->first + 1;
    185:       if (numEntries != PETSC_DETERMINE) {
    186:         /* Allocate space only if we already know the size of the block */
    187:         PetscCall(PetscMalloc1(numEntries * numFaces, (PetscInt **)&s->data));
    188:       }
    189:       for (f = 0; f < numFaces; f++) {
    190:         if (s->nd == 0) {
    191:           /* Determine the size of the block for "mixed" facets */
    192:           PetscInt numFaceVert = 0;
    193:           PetscCall(DMPlexCreateFluent_ReadValues(viewer, &numFaceVert, 1, PETSC_INT, s->index == 2013 ? PETSC_TRUE : PETSC_FALSE));
    194:           if (numEntries == PETSC_DETERMINE) {
    195:             numEntries = numFaceVert + 2;
    196:             PetscCall(PetscMalloc1(numEntries * numFaces, (PetscInt **)&s->data));
    197:           } else {
    198:             PetscCheck(numEntries == numFaceVert + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for mixed faces in Fluent files");
    199:           }
    200:         }
    201:         PetscCall(DMPlexCreateFluent_ReadValues(viewer, &(((PetscInt *)s->data)[f * numEntries]), numEntries, PETSC_INT, s->index == 2013 ? PETSC_TRUE : PETSC_FALSE));
    202:       }
    203:       s->nd = (int)(numEntries - 2);
    204:       PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));
    205:     }
    206:     PetscCall(DMPlexCreateFluent_ReadString(viewer, buffer, ')'));

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

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

    225:   Collective

    227:   Input Parameters:
    228: + comm        - The MPI communicator
    229: . viewer      - The `PetscViewer` associated with a Fluent mesh file
    230: - interpolate - Create faces and edges in the mesh

    232:   Output Parameter:
    233: . dm - The `DM` object representing the mesh

    235:   Level: beginner

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

    251:   PetscFunctionBegin;
    252:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

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

    261:       } else if (s.index == 10 || s.index == 2010) { /* Vertices */
    262:         if (s.zoneID == 0) numVertices = s.last;
    263:         else {
    264:           PetscCheck(!coordsIn, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files");
    265:           coordsIn = (PetscScalar *)s.data;
    266:         }

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

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

    306:           PetscCheck(numFaceVertices == PETSC_DETERMINE || s.nd == numFaceVertices, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed facets in Fluent files are not supported");
    307:           PetscCheck(numFaces >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file");
    308:           if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd;
    309:           numFaceEntries = numFaceVertices + 2;
    310:           if (!faces) PetscCall(PetscMalloc1(numFaces * numFaceEntries, &faces));
    311:           if (!faceZoneIDs) PetscCall(PetscMalloc1(numFaces, &faceZoneIDs));
    312:           PetscCall(PetscMemcpy(&faces[(s.first - 1) * numFaceEntries], s.data, (s.last - s.first + 1) * numFaceEntries * sizeof(PetscInt)));
    313:           /* Record the zoneID for each face set */
    314:           for (z = s.first - 1; z < s.last; z++) faceZoneIDs[z] = s.zoneID;
    315:           PetscCall(PetscFree(s.data));
    316:         }
    317:       }
    318:     } while (s.index >= 0);
    319:   }
    320:   PetscCallMPI(MPI_Bcast(&dim, 1, MPIU_INT, 0, comm));
    321:   PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension");

    323:   /* Allocate cell-vertex mesh */
    324:   PetscCall(DMCreate(comm, dm));
    325:   PetscCall(DMSetType(*dm, DMPLEX));
    326:   PetscCall(DMSetDimension(*dm, dim));
    327:   PetscCall(DMPlexSetChart(*dm, 0, numCells + numVertices));
    328:   if (rank == 0) {
    329:     PetscCheck(numCells >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown number of cells in Fluent file");
    330:     /* If no cell type was given we assume simplices */
    331:     if (numCellVertices == PETSC_DETERMINE) numCellVertices = numFaceVertices + 1;
    332:     for (c = 0; c < numCells; ++c) PetscCall(DMPlexSetConeSize(*dm, c, numCellVertices));
    333:   }
    334:   PetscCall(DMSetUp(*dm));

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

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

    382:     PetscCall(DMPlexInterpolate(*dm, &idm));
    383:     PetscCall(DMDestroy(dm));
    384:     *dm = idm;
    385:   }

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

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

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

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

    450:   if (rank == 0) {
    451:     PetscCall(PetscFree(cellVertices));
    452:     PetscCall(PetscFree(faces));
    453:     PetscCall(PetscFree(faceZoneIDs));
    454:     PetscCall(PetscFree(coordsIn));
    455:   }
    456:   PetscFunctionReturn(PETSC_SUCCESS);
    457: }