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 @*/ 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 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 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 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 @*/ 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 DMLabel faceSets = NULL; 235 PetscInt d, coordSize; 236 PetscScalar *coords, *coordsIn = NULL; 237 PetscSection coordSection; 238 Vec coordinates; 239 PetscErrorCode ierr; 240 241 PetscFunctionBegin; 242 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 243 244 if (!rank) { 245 FluentSection s; 246 do { 247 ierr = DMPlexCreateFluent_ReadSection(viewer, &s);CHKERRQ(ierr); 248 if (s.index == 2) { /* Dimension */ 249 dim = s.nd; 250 251 } else if (s.index == 10 || s.index == 2010) { /* Vertices */ 252 if (s.zoneID == 0) numVertices = s.last; 253 else { 254 if (coordsIn) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files"); 255 coordsIn = (PetscScalar *) s.data; 256 } 257 258 } else if (s.index == 12 || s.index == 2012) { /* Cells */ 259 if (s.zoneID == 0) numCells = s.last; 260 else { 261 switch (s.nd) { 262 case 0: numCellVertices = PETSC_DETERMINE; break; 263 case 1: numCellVertices = 3; break; /* triangular */ 264 case 2: numCellVertices = 4; break; /* tetrahedral */ 265 case 3: numCellVertices = 4; break; /* quadrilateral */ 266 case 4: numCellVertices = 8; break; /* hexahedral */ 267 case 5: numCellVertices = 5; break; /* pyramid */ 268 case 6: numCellVertices = 6; break; /* wedge */ 269 default: numCellVertices = PETSC_DETERMINE; 270 } 271 } 272 273 } else if (s.index == 13 || s.index == 2013) { /* Facets */ 274 if (s.zoneID == 0) { /* Header section */ 275 numFaces = (PetscInt) (s.last - s.first + 1); 276 if (s.nd == 0 || s.nd == 5) numFaceVertices = PETSC_DETERMINE; 277 else numFaceVertices = s.nd; 278 } else { /* Data section */ 279 unsigned int z; 280 281 if (numFaceVertices != PETSC_DETERMINE && s.nd != numFaceVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed facets in Fluent files are not supported"); 282 if (numFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file"); 283 if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd; 284 numFaceEntries = numFaceVertices + 2; 285 if (!faces) {ierr = PetscMalloc1(numFaces*numFaceEntries, &faces);CHKERRQ(ierr);} 286 if (!faceZoneIDs) {ierr = PetscMalloc1(numFaces, &faceZoneIDs);CHKERRQ(ierr);} 287 ierr = PetscMemcpy(&faces[(s.first-1)*numFaceEntries], s.data, (s.last-s.first+1)*numFaceEntries*sizeof(PetscInt));CHKERRQ(ierr); 288 /* Record the zoneID for each face set */ 289 for (z = s.first -1; z < s.last; z++) faceZoneIDs[z] = s.zoneID; 290 ierr = PetscFree(s.data);CHKERRQ(ierr); 291 } 292 } 293 } while (s.index >= 0); 294 } 295 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRMPI(ierr); 296 if (dim < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension"); 297 298 /* Allocate cell-vertex mesh */ 299 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 300 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 301 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 302 ierr = DMPlexSetChart(*dm, 0, numCells + numVertices);CHKERRQ(ierr); 303 if (!rank) { 304 if (numCells < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown number of cells in Fluent file"); 305 /* If no cell type was given we assume simplices */ 306 if (numCellVertices == PETSC_DETERMINE) numCellVertices = numFaceVertices + 1; 307 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(*dm, c, numCellVertices);CHKERRQ(ierr);} 308 } 309 ierr = DMSetUp(*dm);CHKERRQ(ierr); 310 311 if (!rank && faces) { 312 /* Derive cell-vertex list from face-vertex and face-cell maps */ 313 ierr = PetscMalloc1(numCells*numCellVertices, &cellVertices);CHKERRQ(ierr); 314 for (c = 0; c < numCells*numCellVertices; c++) cellVertices[c] = -1; 315 for (f = 0; f < numFaces; f++) { 316 PetscInt *cell; 317 const PetscInt cl = faces[f*numFaceEntries + numFaceVertices]; 318 const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1]; 319 const PetscInt *face = &(faces[f*numFaceEntries]); 320 321 if (cl > 0) { 322 cell = &(cellVertices[(cl-1) * numCellVertices]); 323 for (v = 0; v < numFaceVertices; v++) { 324 PetscBool found = PETSC_FALSE; 325 for (c = 0; c < numCellVertices; c++) { 326 if (cell[c] < 0) break; 327 if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;} 328 } 329 if (!found) cell[c] = face[v]-1 + numCells; 330 } 331 } 332 if (cr > 0) { 333 cell = &(cellVertices[(cr-1) * numCellVertices]); 334 for (v = 0; v < numFaceVertices; v++) { 335 PetscBool found = PETSC_FALSE; 336 for (c = 0; c < numCellVertices; c++) { 337 if (cell[c] < 0) break; 338 if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;} 339 } 340 if (!found) cell[c] = face[v]-1 + numCells; 341 } 342 } 343 } 344 for (c = 0; c < numCells; c++) { 345 ierr = DMPlexSetCone(*dm, c, &(cellVertices[c*numCellVertices]));CHKERRQ(ierr); 346 } 347 } 348 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 349 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 350 if (interpolate) { 351 DM idm; 352 353 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 354 ierr = DMDestroy(dm);CHKERRQ(ierr); 355 *dm = idm; 356 } 357 358 if (!rank && faces) { 359 PetscInt fi, joinSize, meetSize, *fverts, cells[2]; 360 const PetscInt *join, *meet; 361 ierr = PetscMalloc1(numFaceVertices, &fverts);CHKERRQ(ierr); 362 /* Mark facets by finding the full join of all adjacent vertices */ 363 for (f = 0; f < numFaces; f++) { 364 const PetscInt cl = faces[f*numFaceEntries + numFaceVertices] - 1; 365 const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1] - 1; 366 if (cl > 0 && cr > 0) { 367 /* If we know both adjoining cells we can use a single-level meet */ 368 cells[0] = cl; cells[1] = cr; 369 ierr = DMPlexGetMeet(*dm, 2, cells, &meetSize, &meet);CHKERRQ(ierr); 370 if (meetSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %d", f); 371 ierr = DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", meet[0], faceZoneIDs[f]);CHKERRQ(ierr); 372 ierr = DMPlexRestoreMeet(*dm, numFaceVertices, fverts, &meetSize, &meet);CHKERRQ(ierr); 373 } else { 374 for (fi = 0; fi < numFaceVertices; fi++) fverts[fi] = faces[f*numFaceEntries + fi] + numCells - 1; 375 ierr = DMPlexGetFullJoin(*dm, numFaceVertices, fverts, &joinSize, &join);CHKERRQ(ierr); 376 if (joinSize != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not determine Plex facet for Fluent face %d", f); 377 ierr = DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], faceZoneIDs[f]);CHKERRQ(ierr); 378 ierr = DMPlexRestoreJoin(*dm, numFaceVertices, fverts, &joinSize, &join);CHKERRQ(ierr); 379 } 380 } 381 ierr = PetscFree(fverts);CHKERRQ(ierr); 382 } 383 384 { /* Create Face Sets label at all processes */ 385 enum {n = 1}; 386 PetscBool flag[n]; 387 388 flag[0] = faceSets ? PETSC_TRUE : PETSC_FALSE; 389 ierr = MPI_Bcast(flag, n, MPIU_BOOL, 0, comm);CHKERRMPI(ierr); 390 if (flag[0]) {ierr = DMCreateLabel(*dm, "Face Sets");CHKERRQ(ierr);} 391 } 392 393 /* Read coordinates */ 394 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 395 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 396 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 397 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 398 for (v = numCells; v < numCells+numVertices; ++v) { 399 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 400 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 401 } 402 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 403 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 404 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr); 405 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 406 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 407 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 408 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 409 if (!rank && coordsIn) { 410 for (v = 0; v < numVertices; ++v) { 411 for (d = 0; d < dim; ++d) { 412 coords[v*dim+d] = coordsIn[v*dim+d]; 413 } 414 } 415 } 416 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 417 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 418 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 419 420 if (!rank) { 421 ierr = PetscFree(cellVertices);CHKERRQ(ierr); 422 ierr = PetscFree(faces);CHKERRQ(ierr); 423 ierr = PetscFree(faceZoneIDs);CHKERRQ(ierr); 424 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 425 } 426 PetscFunctionReturn(0); 427 } 428