xref: /petsc/src/dm/impls/plex/plexfluent.c (revision b48ab24ab3d10fe53e9be11bf7e4243edd5ad7ec)
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 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 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 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, numFaceVert;
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           ierr = DMPlexCreateFluent_ReadValues(viewer, &numFaceVert, 1, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);CHKERRQ(ierr);
179           if (numEntries == PETSC_DETERMINE) {
180             numEntries = numFaceVert + 2;
181             ierr = PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);CHKERRQ(ierr);
182           } else {
183             if (numEntries != numFaceVert + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for mixed faces in Fluent files");
184           }
185         }
186         ierr = DMPlexCreateFluent_ReadValues(viewer, &(((PetscInt*)s->data)[f*numEntries]), numEntries, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);CHKERRQ(ierr);
187       }
188       s->nd = numEntries - 2;
189       ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
190     }
191     ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
192 
193   } else {                       /* Unknown section type */
194     PetscInt depth = 1;
195     do {
196       /* Match parentheses when parsing unknown sections */
197       do {ierr = PetscViewerRead(viewer, &(buffer[0]), 1, NULL, PETSC_CHAR);CHKERRQ(ierr);}
198       while (buffer[0] != '(' && buffer[0] != ')');
199       if (buffer[0] == '(') depth++;
200       if (buffer[0] == ')') depth--;
201     } while (depth > 0);
202     ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '\n');CHKERRQ(ierr);
203   }
204   PetscFunctionReturn(0);
205 }
206 
207 /*@C
208   DMPlexCreateFluent - Create a DMPlex mesh from a Fluent mesh file.
209 
210   Collective on comm
211 
212   Input Parameters:
213 + comm  - The MPI communicator
214 . viewer - The Viewer associated with a Fluent mesh file
215 - interpolate - Create faces and edges in the mesh
216 
217   Output Parameter:
218 . dm  - The DM object representing the mesh
219 
220   Note: http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm
221 
222   Level: beginner
223 
224 .keywords: mesh, fluent, case
225 .seealso: DMPLEX, DMCreate()
226 @*/
227 PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
228 {
229   PetscMPIInt    rank;
230   PetscInt       c, v, dim = PETSC_DETERMINE, numCells = 0, numVertices = 0, numCellVertices = PETSC_DETERMINE;
231   PetscInt       numFaces = PETSC_DETERMINE, f, numFaceEntries = PETSC_DETERMINE, numFaceVertices = PETSC_DETERMINE;
232   PetscInt      *faces = NULL, *cellVertices = NULL, *faceZoneIDs = NULL;
233   PetscInt       d, coordSize;
234   PetscScalar   *coords, *coordsIn = NULL;
235   PetscSection   coordSection;
236   Vec            coordinates;
237   PetscErrorCode ierr;
238 
239   PetscFunctionBegin;
240   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
241 
242   if (!rank) {
243     FluentSection s;
244     do {
245       ierr = DMPlexCreateFluent_ReadSection(viewer, &s);CHKERRQ(ierr);
246       if (s.index == 2) {                 /* Dimension */
247         dim = s.nd;
248 
249       } else if (s.index == 10 || s.index == 2010) { /* Vertices */
250         if (s.zoneID == 0) numVertices = s.last;
251         else {
252           if (coordsIn) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files");
253           coordsIn = (PetscScalar *) s.data;
254         }
255 
256       } else if (s.index == 12 || s.index == 2012) { /* Cells */
257         if (s.zoneID == 0) numCells = s.last;
258         else {
259           switch (s.nd) {
260           case 0: numCellVertices = PETSC_DETERMINE; break;
261           case 1: numCellVertices = 3; break;  /* triangular */
262           case 2: numCellVertices = 4; break;  /* tetrahedral */
263           case 3: numCellVertices = 4; break;  /* quadrilateral */
264           case 4: numCellVertices = 8; break;  /* hexahedral */
265           case 5: numCellVertices = 5; break;  /* pyramid */
266           case 6: numCellVertices = 6; break;  /* wedge */
267           default: numCellVertices = PETSC_DETERMINE;
268           }
269         }
270 
271       } else if (s.index == 13 || s.index == 2013) { /* Facets */
272         if (s.zoneID == 0) {  /* Header section */
273           numFaces = (PetscInt) (s.last - s.first + 1);
274           if (s.nd == 0 || s.nd == 5) numFaceVertices = PETSC_DETERMINE;
275           else numFaceVertices = s.nd;
276         } else {              /* Data section */
277           unsigned int z;
278 
279           if (numFaceVertices != PETSC_DETERMINE && s.nd != numFaceVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed facets in Fluent files are not supported");
280           if (numFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file");
281           if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd;
282           numFaceEntries = numFaceVertices + 2;
283           if (!faces) {ierr = PetscMalloc1(numFaces*numFaceEntries, &faces);CHKERRQ(ierr);}
284           if (!faceZoneIDs) {ierr = PetscMalloc1(numFaces, &faceZoneIDs);CHKERRQ(ierr);}
285           ierr = PetscMemcpy(&(faces[(s.first-1)*numFaceEntries]), s.data, (s.last-s.first+1)*numFaceEntries*sizeof(PetscInt));CHKERRQ(ierr);
286           /* Record the zoneID for each face set */
287           for (z = s.first -1; z < s.last; z++) faceZoneIDs[z] = s.zoneID;
288           ierr = PetscFree(s.data);CHKERRQ(ierr);
289         }
290       }
291     } while (s.index >= 0);
292   }
293   ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr);
294   if (dim < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension");
295 
296   /* Allocate cell-vertex mesh */
297   ierr = DMCreate(comm, dm);CHKERRQ(ierr);
298   ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
299   ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
300   ierr = DMPlexSetChart(*dm, 0, numCells + numVertices);CHKERRQ(ierr);
301   if (!rank) {
302     if (numCells < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown number of cells in Fluent file");
303     /* If no cell type was given we assume simplices */
304     if (numCellVertices == PETSC_DETERMINE) numCellVertices = numFaceVertices + 1;
305     for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(*dm, c, numCellVertices);CHKERRQ(ierr);}
306   }
307   ierr = DMSetUp(*dm);CHKERRQ(ierr);
308 
309   if (!rank && faces) {
310     /* Derive cell-vertex list from face-vertex and face-cell maps */
311     ierr = PetscMalloc1(numCells*numCellVertices, &cellVertices);CHKERRQ(ierr);
312     for (c = 0; c < numCells*numCellVertices; c++) cellVertices[c] = -1;
313     for (f = 0; f < numFaces; f++) {
314       PetscInt *cell;
315       const PetscInt cl = faces[f*numFaceEntries + numFaceVertices];
316       const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1];
317       const PetscInt *face = &(faces[f*numFaceEntries]);
318 
319       if (cl > 0) {
320         cell = &(cellVertices[(cl-1) * numCellVertices]);
321         for (v = 0; v < numFaceVertices; v++) {
322           PetscBool found = PETSC_FALSE;
323           for (c = 0; c < numCellVertices; c++) {
324             if (cell[c] < 0) break;
325             if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
326           }
327           if (!found) cell[c] = face[v]-1 + numCells;
328         }
329       }
330       if (cr > 0) {
331         cell = &(cellVertices[(cr-1) * numCellVertices]);
332         for (v = 0; v < numFaceVertices; v++) {
333           PetscBool found = PETSC_FALSE;
334           for (c = 0; c < numCellVertices; c++) {
335             if (cell[c] < 0) break;
336             if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;}
337           }
338           if (!found) cell[c] = face[v]-1 + numCells;
339         }
340       }
341     }
342     for (c = 0; c < numCells; c++) {
343       ierr = DMPlexSetCone(*dm, c, &(cellVertices[c*numCellVertices]));CHKERRQ(ierr);
344     }
345   }
346   ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
347   ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
348   if (interpolate) {
349     DM idm = NULL;
350 
351     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
352     ierr = DMDestroy(dm);CHKERRQ(ierr);
353     *dm  = idm;
354   }
355 
356   if (!rank && faces) {
357     PetscInt fi, joinSize, meetSize, *fverts, cells[2];
358     const PetscInt *join, *meet;
359     ierr = PetscMalloc1(numFaceVertices, &fverts);CHKERRQ(ierr);
360     /* Mark facets by finding the full join of all adjacent vertices */
361     for (f = 0; f < numFaces; f++) {
362       const PetscInt cl = faces[f*numFaceEntries + numFaceVertices] - 1;
363       const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1] - 1;
364       if (cl > 0 && cr > 0) {
365         /* If we know both adjoining cells we can use a single-level meet */
366         cells[0] = cl; cells[1] = cr;
367         ierr = DMPlexGetMeet(*dm, 2, cells, &meetSize, &meet);CHKERRQ(ierr);
368         if (meetSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %d", f);
369         ierr = DMSetLabelValue(*dm, "Face Sets", meet[0], faceZoneIDs[f]);CHKERRQ(ierr);
370         ierr = DMPlexRestoreMeet(*dm, numFaceVertices, fverts, &meetSize, &meet);CHKERRQ(ierr);
371       } else {
372         for (fi = 0; fi < numFaceVertices; fi++) fverts[fi] = faces[f*numFaceEntries + fi] + numCells - 1;
373         ierr = DMPlexGetFullJoin(*dm, numFaceVertices, fverts, &joinSize, &join);CHKERRQ(ierr);
374         if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %d", f);
375         ierr = DMSetLabelValue(*dm, "Face Sets", join[0], faceZoneIDs[f]);CHKERRQ(ierr);
376         ierr = DMPlexRestoreJoin(*dm, numFaceVertices, fverts, &joinSize, &join);CHKERRQ(ierr);
377       }
378     }
379     ierr = PetscFree(fverts);CHKERRQ(ierr);
380   }
381 
382   /* Read coordinates */
383   ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
384   ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
385   ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr);
386   ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr);
387   for (v = numCells; v < numCells+numVertices; ++v) {
388     ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr);
389     ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr);
390   }
391   ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
392   ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
393   ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
394   ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
395   ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
396   ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr);
397   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
398   if (!rank && coordsIn) {
399     for (v = 0; v < numVertices; ++v) {
400       for (d = 0; d < dim; ++d) {
401         coords[v*dim+d] = coordsIn[v*dim+d];
402       }
403     }
404   }
405   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
406   ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr);
407   ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
408   if (!rank) {
409     if (cellVertices) {ierr = PetscFree(cellVertices);CHKERRQ(ierr);}
410     ierr = PetscFree(faces);CHKERRQ(ierr);
411     ierr = PetscFree(faceZoneIDs);CHKERRQ(ierr);
412     ierr = PetscFree(coordsIn);CHKERRQ(ierr);
413   }
414   PetscFunctionReturn(0);
415 }
416