xref: /petsc/src/dm/impls/plex/plexfluent.c (revision b4db17ed4967172e430e8d3533aad9efd6fc38f9)
1 #define PETSCDM_DLL
2 #include <petsc/private/dmpleximpl.h>    /*I   "petscdmplex.h"   I*/
3 
4 /*@C
5   DMPlexCreateFluentFromFile - Create a DMPlex mesh from a Fluent mesh file
6 
7 + comm        - The MPI communicator
8 . filename    - Name of the Fluent mesh file
9 - interpolate - Create faces and edges in the mesh
10 
11   Output Parameter:
12 . dm  - The DM object representing the mesh
13 
14   Level: beginner
15 
16 .seealso: DMPlexCreateFromFile(), DMPlexCreateFluent(), DMPlexCreate()
17 @*/
18 PetscErrorCode DMPlexCreateFluentFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
19 {
20   PetscViewer     viewer;
21   PetscErrorCode  ierr;
22 
23   PetscFunctionBegin;
24   /* Create file viewer and build plex */
25   ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr);
26   ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr);
27   ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr);
28   ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr);
29   ierr = DMPlexCreateFluent(comm, viewer, interpolate, dm);CHKERRQ(ierr);
30   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
31   PetscFunctionReturn(0);
32 }
33 
34 static PetscErrorCode DMPlexCreateFluent_ReadString(PetscViewer viewer, char *buffer, char delim)
35 {
36   PetscInt ret, i = 0;
37   PetscErrorCode ierr;
38 
39   PetscFunctionBegin;
40   do {ierr = PetscViewerRead(viewer, &(buffer[i++]), 1, &ret, PETSC_CHAR);CHKERRQ(ierr);}
41   while (ret > 0 && buffer[i-1] != '\0' && buffer[i-1] != delim);
42   buffer[i] = '\0';
43   PetscFunctionReturn(0);
44 }
45 
46 static PetscErrorCode DMPlexCreateFluent_ReadValues(PetscViewer viewer, void *data, PetscInt count, PetscDataType dtype, PetscBool binary)
47 {
48   int            fdes=0;
49   FILE          *file;
50   PetscInt       i;
51   PetscErrorCode ierr;
52 
53   PetscFunctionBegin;
54   if (binary) {
55     /* Extract raw file descriptor to read binary block */
56     ierr = PetscViewerASCIIGetPointer(viewer, &file);CHKERRQ(ierr);
57     fflush(file); fdes = fileno(file);
58   }
59 
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       ierr = PetscViewerRead(viewer, cbuf, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
67       snum = sscanf(cbuf, "%x", &ibuf);
68       if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
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     ierr = PetscMalloc1(count, &ibuf);CHKERRQ(ierr);
75     ierr = PetscBinaryRead(fdes, ibuf, count, PETSC_ENUM);CHKERRQ(ierr);
76     ierr = PetscByteSwap(ibuf, PETSC_ENUM, count);CHKERRQ(ierr);
77     for (i = 0; i < count; i++) ((PetscInt*)data)[i] = (PetscInt)(ibuf[i]);
78     ierr = PetscFree(ibuf);CHKERRQ(ierr);
79 
80  } else if (binary && dtype == PETSC_SCALAR) {
81     float *fbuf;
82     /* Always read 32-bit floats and cast to PetscScalar */
83     ierr = PetscMalloc1(count, &fbuf);CHKERRQ(ierr);
84     ierr = PetscBinaryRead(fdes, fbuf, count, PETSC_FLOAT);CHKERRQ(ierr);
85     ierr = PetscByteSwap(fbuf, PETSC_FLOAT, count);CHKERRQ(ierr);
86     for (i = 0; i < count; i++) ((PetscScalar*)data)[i] = (PetscScalar)(fbuf[i]);
87     ierr = PetscFree(fbuf);CHKERRQ(ierr);
88   } else {
89     ierr = PetscViewerASCIIRead(viewer, data, count, NULL, dtype);CHKERRQ(ierr);
90   }
91   PetscFunctionReturn(0);
92 }
93 
94 static PetscErrorCode DMPlexCreateFluent_ReadSection(PetscViewer viewer, FluentSection *s)
95 {
96   char           buffer[PETSC_MAX_PATH_LEN];
97   int            snum;
98   PetscErrorCode ierr;
99 
100   PetscFunctionBegin;
101   /* Fast-forward to next section and derive its index */
102   ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr);
103   ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ' ');CHKERRQ(ierr);
104   snum = sscanf(buffer, "%d", &(s->index));
105   /* If we can't match an index return -1 to signal end-of-file */
106   if (snum < 1) {s->index = -1;   PetscFunctionReturn(0);}
107 
108   if (s->index == 0) {           /* Comment */
109     ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
110 
111   } else if (s->index == 2) {    /* Dimension */
112     ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
113     snum = sscanf(buffer, "%d", &(s->nd));
114     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
115 
116   } else if (s->index == 10 || s->index == 2010) {   /* Vertices */
117     ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
118     snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
119     if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
120     if (s->zoneID > 0) {
121       PetscInt numCoords = s->last - s->first + 1;
122       ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr);
123       ierr = PetscMalloc1(s->nd*numCoords, (PetscScalar**)&s->data);CHKERRQ(ierr);
124       ierr = DMPlexCreateFluent_ReadValues(viewer, s->data, s->nd*numCoords, PETSC_SCALAR, s->index==2010 ? PETSC_TRUE : PETSC_FALSE);CHKERRQ(ierr);
125       ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
126     }
127     ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
128 
129   } else if (s->index == 12 || s->index == 2012) {   /* Cells */
130     ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
131     snum = sscanf(buffer, "(%x", &(s->zoneID));
132     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
133     if (s->zoneID == 0) {  /* Header section */
134       snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
135       if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
136     } else {               /* Data section */
137       snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
138       if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
139       if (s->nd == 0) {
140         /* Read cell type definitions for mixed cells */
141         PetscInt numCells = s->last - s->first + 1;
142         ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr);
143         ierr = PetscMalloc1(numCells, (PetscInt**)&s->data);CHKERRQ(ierr);
144         ierr = DMPlexCreateFluent_ReadValues(viewer, s->data, numCells, PETSC_INT, s->index==2012 ? PETSC_TRUE : PETSC_FALSE);CHKERRQ(ierr);
145         ierr = PetscFree(s->data);CHKERRQ(ierr);
146         ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
147       }
148     }
149     ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
150 
151   } else if (s->index == 13 || s->index == 2013) {   /* Faces */
152     ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
153     snum = sscanf(buffer, "(%x", &(s->zoneID));
154     if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
155     if (s->zoneID == 0) {  /* Header section */
156       snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
157       if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
158     } else {               /* Data section */
159       PetscInt f, numEntries, numFaces;
160       snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
161       if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
162       ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr);
163       switch (s->nd) {
164       case 0: numEntries = PETSC_DETERMINE; break;
165       case 2: numEntries = 2 + 2; break;  /* linear */
166       case 3: numEntries = 2 + 3; break;  /* triangular */
167       case 4: numEntries = 2 + 4; break;  /* quadrilateral */
168       default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file");
169       }
170       numFaces = s->last-s->first + 1;
171       if (numEntries != PETSC_DETERMINE) {
172         /* Allocate space only if we already know the size of the block */
173         ierr = PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);CHKERRQ(ierr);
174       }
175       for (f = 0; f < numFaces; f++) {
176         if (s->nd == 0) {
177           /* Determine the size of the block for "mixed" facets */
178           PetscInt numFaceVert = 0;
179           ierr = DMPlexCreateFluent_ReadValues(viewer, &numFaceVert, 1, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);CHKERRQ(ierr);
180           if (numEntries == PETSC_DETERMINE) {
181             numEntries = numFaceVert + 2;
182             ierr = PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);CHKERRQ(ierr);
183           } else {
184             if (numEntries != numFaceVert + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for mixed faces in Fluent files");
185           }
186         }
187         ierr = DMPlexCreateFluent_ReadValues(viewer, &(((PetscInt*)s->data)[f*numEntries]), numEntries, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);CHKERRQ(ierr);
188       }
189       s->nd = numEntries - 2;
190       ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
191     }
192     ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
193 
194   } else {                       /* Unknown section type */
195     PetscInt depth = 1;
196     do {
197       /* Match parentheses when parsing unknown sections */
198       do {ierr = PetscViewerRead(viewer, &(buffer[0]), 1, NULL, PETSC_CHAR);CHKERRQ(ierr);}
199       while (buffer[0] != '(' && buffer[0] != ')');
200       if (buffer[0] == '(') depth++;
201       if (buffer[0] == ')') depth--;
202     } while (depth > 0);
203     ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '\n');CHKERRQ(ierr);
204   }
205   PetscFunctionReturn(0);
206 }
207 
208 /*@C
209   DMPlexCreateFluent - Create a DMPlex mesh from a Fluent mesh file.
210 
211   Collective on comm
212 
213   Input Parameters:
214 + comm  - The MPI communicator
215 . viewer - The Viewer associated with a Fluent mesh file
216 - interpolate - Create faces and edges in the mesh
217 
218   Output Parameter:
219 . dm  - The DM object representing the mesh
220 
221   Note: http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm
222 
223   Level: beginner
224 
225 .keywords: mesh, fluent, case
226 .seealso: DMPLEX, DMCreate()
227 @*/
228 PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
229 {
230   PetscMPIInt    rank;
231   PetscInt       c, v, dim = PETSC_DETERMINE, numCells = 0, numVertices = 0, numCellVertices = PETSC_DETERMINE;
232   PetscInt       numFaces = PETSC_DETERMINE, f, numFaceEntries = PETSC_DETERMINE, numFaceVertices = PETSC_DETERMINE;
233   PetscInt      *faces = NULL, *cellVertices = NULL, *faceZoneIDs = NULL;
234   PetscInt       d, coordSize;
235   PetscScalar   *coords, *coordsIn = NULL;
236   PetscSection   coordSection;
237   Vec            coordinates;
238   PetscErrorCode ierr;
239 
240   PetscFunctionBegin;
241   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
242 
243   if (!rank) {
244     FluentSection s;
245     do {
246       ierr = DMPlexCreateFluent_ReadSection(viewer, &s);CHKERRQ(ierr);
247       if (s.index == 2) {                 /* Dimension */
248         dim = s.nd;
249 
250       } else if (s.index == 10 || s.index == 2010) { /* Vertices */
251         if (s.zoneID == 0) numVertices = s.last;
252         else {
253           if (coordsIn) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files");
254           coordsIn = (PetscScalar *) s.data;
255         }
256 
257       } else if (s.index == 12 || s.index == 2012) { /* Cells */
258         if (s.zoneID == 0) numCells = s.last;
259         else {
260           switch (s.nd) {
261           case 0: numCellVertices = PETSC_DETERMINE; break;
262           case 1: numCellVertices = 3; break;  /* triangular */
263           case 2: numCellVertices = 4; break;  /* tetrahedral */
264           case 3: numCellVertices = 4; break;  /* quadrilateral */
265           case 4: numCellVertices = 8; break;  /* hexahedral */
266           case 5: numCellVertices = 5; break;  /* pyramid */
267           case 6: numCellVertices = 6; break;  /* wedge */
268           default: numCellVertices = PETSC_DETERMINE;
269           }
270         }
271 
272       } else if (s.index == 13 || s.index == 2013) { /* Facets */
273         if (s.zoneID == 0) {  /* Header section */
274           numFaces = (PetscInt) (s.last - s.first + 1);
275           if (s.nd == 0 || s.nd == 5) numFaceVertices = PETSC_DETERMINE;
276           else numFaceVertices = s.nd;
277         } else {              /* Data section */
278           unsigned int z;
279 
280           if (numFaceVertices != PETSC_DETERMINE && s.nd != numFaceVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed facets in Fluent files are not supported");
281           if (numFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file");
282           if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd;
283           numFaceEntries = numFaceVertices + 2;
284           if (!faces) {ierr = PetscMalloc1(numFaces*numFaceEntries, &faces);CHKERRQ(ierr);}
285           if (!faceZoneIDs) {ierr = PetscMalloc1(numFaces, &faceZoneIDs);CHKERRQ(ierr);}
286           ierr = PetscMemcpy(&(faces[(s.first-1)*numFaceEntries]), s.data, (s.last-s.first+1)*numFaceEntries*sizeof(PetscInt));CHKERRQ(ierr);
287           /* Record the zoneID for each face set */
288           for (z = s.first -1; z < s.last; z++) faceZoneIDs[z] = s.zoneID;
289           ierr = PetscFree(s.data);CHKERRQ(ierr);
290         }
291       }
292     } while (s.index >= 0);
293   }
294   ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr);
295   if (dim < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension");
296 
297   /* Allocate cell-vertex mesh */
298   ierr = DMCreate(comm, dm);CHKERRQ(ierr);
299   ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
300   ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
301   ierr = DMPlexSetChart(*dm, 0, numCells + numVertices);CHKERRQ(ierr);
302   if (!rank) {
303     if (numCells < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown number of cells in Fluent file");
304     /* If no cell type was given we assume simplices */
305     if (numCellVertices == PETSC_DETERMINE) numCellVertices = numFaceVertices + 1;
306     for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(*dm, c, numCellVertices);CHKERRQ(ierr);}
307   }
308   ierr = DMSetUp(*dm);CHKERRQ(ierr);
309 
310   if (!rank && faces) {
311     /* Derive cell-vertex list from face-vertex and face-cell maps */
312     ierr = PetscMalloc1(numCells*numCellVertices, &cellVertices);CHKERRQ(ierr);
313     for (c = 0; c < numCells*numCellVertices; c++) cellVertices[c] = -1;
314     for (f = 0; f < numFaces; f++) {
315       PetscInt *cell;
316       const PetscInt cl = faces[f*numFaceEntries + numFaceVertices];
317       const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1];
318       const PetscInt *face = &(faces[f*numFaceEntries]);
319 
320       if (cl > 0) {
321         cell = &(cellVertices[(cl-1) * numCellVertices]);
322         for (v = 0; v < numFaceVertices; v++) {
323           PetscBool found = PETSC_FALSE;
324           for (c = 0; c < numCellVertices; c++) {
325             if (cell[c] < 0) break;
326             if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
327           }
328           if (!found) cell[c] = face[v]-1 + numCells;
329         }
330       }
331       if (cr > 0) {
332         cell = &(cellVertices[(cr-1) * numCellVertices]);
333         for (v = 0; v < numFaceVertices; v++) {
334           PetscBool found = PETSC_FALSE;
335           for (c = 0; c < numCellVertices; c++) {
336             if (cell[c] < 0) break;
337             if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
338           }
339           if (!found) cell[c] = face[v]-1 + numCells;
340         }
341       }
342     }
343     for (c = 0; c < numCells; c++) {
344       ierr = DMPlexSetCone(*dm, c, &(cellVertices[c*numCellVertices]));CHKERRQ(ierr);
345     }
346   }
347   ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
348   ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
349   if (interpolate) {
350     DM idm;
351 
352     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
353     ierr = DMDestroy(dm);CHKERRQ(ierr);
354     *dm  = idm;
355   }
356 
357   if (!rank && faces) {
358     PetscInt fi, joinSize, meetSize, *fverts, cells[2];
359     const PetscInt *join, *meet;
360     ierr = PetscMalloc1(numFaceVertices, &fverts);CHKERRQ(ierr);
361     /* Mark facets by finding the full join of all adjacent vertices */
362     for (f = 0; f < numFaces; f++) {
363       const PetscInt cl = faces[f*numFaceEntries + numFaceVertices] - 1;
364       const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1] - 1;
365       if (cl > 0 && cr > 0) {
366         /* If we know both adjoining cells we can use a single-level meet */
367         cells[0] = cl; cells[1] = cr;
368         ierr = DMPlexGetMeet(*dm, 2, cells, &meetSize, &meet);CHKERRQ(ierr);
369         if (meetSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %d", f);
370         ierr = DMSetLabelValue(*dm, "Face Sets", meet[0], faceZoneIDs[f]);CHKERRQ(ierr);
371         ierr = DMPlexRestoreMeet(*dm, numFaceVertices, fverts, &meetSize, &meet);CHKERRQ(ierr);
372       } else {
373         for (fi = 0; fi < numFaceVertices; fi++) fverts[fi] = faces[f*numFaceEntries + fi] + numCells - 1;
374         ierr = DMPlexGetFullJoin(*dm, numFaceVertices, fverts, &joinSize, &join);CHKERRQ(ierr);
375         if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %d", f);
376         ierr = DMSetLabelValue(*dm, "Face Sets", join[0], faceZoneIDs[f]);CHKERRQ(ierr);
377         ierr = DMPlexRestoreJoin(*dm, numFaceVertices, fverts, &joinSize, &join);CHKERRQ(ierr);
378       }
379     }
380     ierr = PetscFree(fverts);CHKERRQ(ierr);
381   }
382 
383   /* Read coordinates */
384   ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
385   ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
386   ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr);
387   ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr);
388   for (v = numCells; v < numCells+numVertices; ++v) {
389     ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr);
390     ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr);
391   }
392   ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
393   ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
394   ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
395   ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
396   ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
397   ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr);
398   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
399   if (!rank && 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   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
407   ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr);
408   ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
409   if (!rank) {
410     if (cellVertices) {ierr = PetscFree(cellVertices);CHKERRQ(ierr);}
411     ierr = PetscFree(faces);CHKERRQ(ierr);
412     ierr = PetscFree(faceZoneIDs);CHKERRQ(ierr);
413     ierr = PetscFree(coordsIn);CHKERRQ(ierr);
414   }
415   PetscFunctionReturn(0);
416 }
417