1 #define PETSC_DESIRE_FEATURE_TEST_MACROS /* for fileno() */
2 #define PETSCDM_DLL
3 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
4
5 /*@C
6 DMPlexCreateFluentFromFile - Create a DMPlex mesh from a Fluent mesh file
7
8 + comm - The MPI communicator
9 . filename - Name of the Fluent mesh file
10 - interpolate - Create faces and edges in the mesh
11
12 Output Parameter:
13 . dm - The DM object representing the mesh
14
15 Level: beginner
16
17 .seealso: DMPlexCreateFromFile(), DMPlexCreateFluent(), DMPlexCreate()
18 @*/
DMPlexCreateFluentFromFile(MPI_Comm comm,const char filename[],PetscBool interpolate,DM * dm)19 PetscErrorCode DMPlexCreateFluentFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
20 {
21 PetscViewer viewer;
22 PetscErrorCode ierr;
23
24 PetscFunctionBegin;
25 /* Create file viewer and build plex */
26 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr);
27 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr);
28 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr);
29 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr);
30 ierr = DMPlexCreateFluent(comm, viewer, interpolate, dm);CHKERRQ(ierr);
31 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
32 PetscFunctionReturn(0);
33 }
34
DMPlexCreateFluent_ReadString(PetscViewer viewer,char * buffer,char delim)35 static PetscErrorCode DMPlexCreateFluent_ReadString(PetscViewer viewer, char *buffer, char delim)
36 {
37 PetscInt ret, i = 0;
38 PetscErrorCode ierr;
39
40 PetscFunctionBegin;
41 do {ierr = PetscViewerRead(viewer, &(buffer[i++]), 1, &ret, PETSC_CHAR);CHKERRQ(ierr);}
42 while (ret > 0 && buffer[i-1] != '\0' && buffer[i-1] != delim);
43 if (!ret) buffer[i-1] = '\0'; else buffer[i] = '\0';
44 PetscFunctionReturn(0);
45 }
46
DMPlexCreateFluent_ReadValues(PetscViewer viewer,void * data,PetscInt count,PetscDataType dtype,PetscBool binary)47 static PetscErrorCode DMPlexCreateFluent_ReadValues(PetscViewer viewer, void *data, PetscInt count, PetscDataType dtype, PetscBool binary)
48 {
49 int fdes=0;
50 FILE *file;
51 PetscInt i;
52 PetscErrorCode ierr;
53
54 PetscFunctionBegin;
55 if (binary) {
56 /* Extract raw file descriptor to read binary block */
57 ierr = PetscViewerASCIIGetPointer(viewer, &file);CHKERRQ(ierr);
58 fflush(file); fdes = fileno(file);
59 }
60
61 if (!binary && dtype == PETSC_INT) {
62 char cbuf[256];
63 unsigned int ibuf;
64 int snum;
65 /* Parse hexadecimal ascii integers */
66 for (i = 0; i < count; i++) {
67 ierr = PetscViewerRead(viewer, cbuf, 1, NULL, PETSC_STRING);CHKERRQ(ierr);
68 snum = sscanf(cbuf, "%x", &ibuf);
69 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
70 ((PetscInt*)data)[i] = (PetscInt)ibuf;
71 }
72 } else if (binary && dtype == PETSC_INT) {
73 /* Always read 32-bit ints and cast to PetscInt */
74 int *ibuf;
75 ierr = PetscMalloc1(count, &ibuf);CHKERRQ(ierr);
76 ierr = PetscBinaryRead(fdes, ibuf, count, NULL, PETSC_ENUM);CHKERRQ(ierr);
77 ierr = PetscByteSwap(ibuf, PETSC_ENUM, count);CHKERRQ(ierr);
78 for (i = 0; i < count; i++) ((PetscInt*)data)[i] = (PetscInt)(ibuf[i]);
79 ierr = PetscFree(ibuf);CHKERRQ(ierr);
80
81 } else if (binary && dtype == PETSC_SCALAR) {
82 float *fbuf;
83 /* Always read 32-bit floats and cast to PetscScalar */
84 ierr = PetscMalloc1(count, &fbuf);CHKERRQ(ierr);
85 ierr = PetscBinaryRead(fdes, fbuf, count, NULL, PETSC_FLOAT);CHKERRQ(ierr);
86 ierr = PetscByteSwap(fbuf, PETSC_FLOAT, count);CHKERRQ(ierr);
87 for (i = 0; i < count; i++) ((PetscScalar*)data)[i] = (PetscScalar)(fbuf[i]);
88 ierr = PetscFree(fbuf);CHKERRQ(ierr);
89 } else {
90 ierr = PetscViewerASCIIRead(viewer, data, count, NULL, dtype);CHKERRQ(ierr);
91 }
92 PetscFunctionReturn(0);
93 }
94
DMPlexCreateFluent_ReadSection(PetscViewer viewer,FluentSection * s)95 static PetscErrorCode DMPlexCreateFluent_ReadSection(PetscViewer viewer, FluentSection *s)
96 {
97 char buffer[PETSC_MAX_PATH_LEN];
98 int snum;
99 PetscErrorCode ierr;
100
101 PetscFunctionBegin;
102 /* Fast-forward to next section and derive its index */
103 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr);
104 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ' ');CHKERRQ(ierr);
105 snum = sscanf(buffer, "%d", &(s->index));
106 /* If we can't match an index return -1 to signal end-of-file */
107 if (snum < 1) {s->index = -1; PetscFunctionReturn(0);}
108
109 if (s->index == 0) { /* Comment */
110 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
111
112 } else if (s->index == 2) { /* Dimension */
113 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
114 snum = sscanf(buffer, "%d", &(s->nd));
115 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
116
117 } else if (s->index == 10 || s->index == 2010) { /* Vertices */
118 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
119 snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
120 if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
121 if (s->zoneID > 0) {
122 PetscInt numCoords = s->last - s->first + 1;
123 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr);
124 ierr = PetscMalloc1(s->nd*numCoords, (PetscScalar**)&s->data);CHKERRQ(ierr);
125 ierr = DMPlexCreateFluent_ReadValues(viewer, s->data, s->nd*numCoords, PETSC_SCALAR, s->index==2010 ? PETSC_TRUE : PETSC_FALSE);CHKERRQ(ierr);
126 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
127 }
128 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
129
130 } else if (s->index == 12 || s->index == 2012) { /* Cells */
131 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
132 snum = sscanf(buffer, "(%x", &(s->zoneID));
133 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
134 if (s->zoneID == 0) { /* Header section */
135 snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
136 if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
137 } else { /* Data section */
138 snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
139 if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
140 if (s->nd == 0) {
141 /* Read cell type definitions for mixed cells */
142 PetscInt numCells = s->last - s->first + 1;
143 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr);
144 ierr = PetscMalloc1(numCells, (PetscInt**)&s->data);CHKERRQ(ierr);
145 ierr = DMPlexCreateFluent_ReadValues(viewer, s->data, numCells, PETSC_INT, s->index==2012 ? PETSC_TRUE : PETSC_FALSE);CHKERRQ(ierr);
146 ierr = PetscFree(s->data);CHKERRQ(ierr);
147 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
148 }
149 }
150 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
151
152 } else if (s->index == 13 || s->index == 2013) { /* Faces */
153 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
154 snum = sscanf(buffer, "(%x", &(s->zoneID));
155 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
156 if (s->zoneID == 0) { /* Header section */
157 snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd));
158 if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
159 } else { /* Data section */
160 PetscInt f, numEntries, numFaces;
161 snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd));
162 if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file");
163 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr);
164 switch (s->nd) {
165 case 0: numEntries = PETSC_DETERMINE; break;
166 case 2: numEntries = 2 + 2; break; /* linear */
167 case 3: numEntries = 2 + 3; break; /* triangular */
168 case 4: numEntries = 2 + 4; break; /* quadrilateral */
169 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file");
170 }
171 numFaces = s->last-s->first + 1;
172 if (numEntries != PETSC_DETERMINE) {
173 /* Allocate space only if we already know the size of the block */
174 ierr = PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);CHKERRQ(ierr);
175 }
176 for (f = 0; f < numFaces; f++) {
177 if (s->nd == 0) {
178 /* Determine the size of the block for "mixed" facets */
179 PetscInt numFaceVert = 0;
180 ierr = DMPlexCreateFluent_ReadValues(viewer, &numFaceVert, 1, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);CHKERRQ(ierr);
181 if (numEntries == PETSC_DETERMINE) {
182 numEntries = numFaceVert + 2;
183 ierr = PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);CHKERRQ(ierr);
184 } else {
185 if (numEntries != numFaceVert + 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for mixed faces in Fluent files");
186 }
187 }
188 ierr = DMPlexCreateFluent_ReadValues(viewer, &(((PetscInt*)s->data)[f*numEntries]), numEntries, PETSC_INT, s->index==2013 ? PETSC_TRUE : PETSC_FALSE);CHKERRQ(ierr);
189 }
190 s->nd = numEntries - 2;
191 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
192 }
193 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr);
194
195 } else { /* Unknown section type */
196 PetscInt depth = 1;
197 do {
198 /* Match parentheses when parsing unknown sections */
199 do {ierr = PetscViewerRead(viewer, &(buffer[0]), 1, NULL, PETSC_CHAR);CHKERRQ(ierr);}
200 while (buffer[0] != '(' && buffer[0] != ')');
201 if (buffer[0] == '(') depth++;
202 if (buffer[0] == ')') depth--;
203 } while (depth > 0);
204 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '\n');CHKERRQ(ierr);
205 }
206 PetscFunctionReturn(0);
207 }
208
209 /*@C
210 DMPlexCreateFluent - Create a DMPlex mesh from a Fluent mesh file.
211
212 Collective
213
214 Input Parameters:
215 + comm - The MPI communicator
216 . viewer - The Viewer associated with a Fluent mesh file
217 - interpolate - Create faces and edges in the mesh
218
219 Output Parameter:
220 . dm - The DM object representing the mesh
221
222 Note: http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm
223
224 Level: beginner
225
226 .seealso: DMPLEX, DMCreate()
227 @*/
DMPlexCreateFluent(MPI_Comm comm,PetscViewer viewer,PetscBool interpolate,DM * dm)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