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