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; 41 PetscErrorCode ierr; 42 43 PetscFunctionBegin; 44 do {ierr = PetscViewerRead(viewer, &(buffer[i++]), 1, 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_ReadSection" 52 PetscErrorCode DMPlexCreateFluent_ReadSection(PetscViewer viewer, FluentSection *s) 53 { 54 char buffer[PETSC_MAX_PATH_LEN]; 55 int snum; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 /* Fast-forward to next section and derive its index */ 60 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr); 61 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ' ');CHKERRQ(ierr); 62 snum = sscanf(buffer, "%d", &(s->index)); 63 /* If we can't match an index return -1 to signal end-of-file */ 64 if (snum < 1) {s->index = -1; PetscFunctionReturn(0);} 65 66 if (s->index == 0) { /* Comment */ 67 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr); 68 69 } else if (s->index == 2) { /* Dimension */ 70 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr); 71 snum = sscanf(buffer, "%d", &(s->nd)); 72 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 73 74 } else if (s->index == 10) { /* Vertices */ 75 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr); 76 snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd)); 77 if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 78 if (s->zoneID > 0) { 79 PetscInt numCoords = s->last - s->first + 1; 80 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr); 81 ierr = PetscMalloc1(s->nd*numCoords, (PetscScalar**)&s->data);CHKERRQ(ierr); 82 ierr = PetscViewerRead(viewer, (PetscScalar*)s->data, numCoords*s->nd, PETSC_REAL);CHKERRQ(ierr); 83 } 84 85 } else if (s->index == 12) { /* Cells */ 86 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr); 87 snum = sscanf(buffer, "(%x", &(s->zoneID)); 88 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 89 if (s->zoneID == 0) { /* Header section */ 90 snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd)); 91 if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 92 } else { /* Data section */ 93 snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd)); 94 if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 95 if (s->nd == 0) { 96 /* Read cell type definitions for mixed cells */ 97 PetscInt i, numCells = s->last - s->first + 1; 98 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr); 99 for (i = 0; i < numCells; i++) { 100 ierr = PetscViewerRead(viewer, buffer, 1, PETSC_STRING);CHKERRQ(ierr); 101 } 102 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr); 103 } 104 } 105 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr); 106 107 } else if (s->index == 13) { /* Faces */ 108 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr); 109 snum = sscanf(buffer, "(%x", &(s->zoneID)); 110 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 111 if (s->zoneID == 0) { /* Header section */ 112 snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd)); 113 if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 114 } else { /* Data section */ 115 PetscInt f, e, numEntries, numFaces; 116 int numFaceVert; 117 snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd)); 118 if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 119 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr); 120 switch (s->nd) { 121 case 0: numEntries = PETSC_DETERMINE; break; 122 case 2: numEntries = 2 + 2; break; /* linear */ 123 case 3: numEntries = 2 + 3; break; /* triangular */ 124 case 4: numEntries = 2 + 4; break; /* quadrilateral */ 125 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file"); 126 } 127 numFaces = s->last-s->first + 1; 128 if (numEntries != PETSC_DETERMINE) { 129 /* Allocate space only if we already know the size of the block */ 130 ierr = PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);CHKERRQ(ierr); 131 } 132 for (f = 0; f < numFaces; f++) { 133 if (s->nd == 0) { 134 /* Determine the size of the block for "mixed" facets */ 135 ierr = PetscViewerRead(viewer, buffer, 1, PETSC_STRING);CHKERRQ(ierr); 136 snum = sscanf(buffer, "%x", &numFaceVert); 137 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 138 if (numEntries == PETSC_DETERMINE) { 139 numEntries = numFaceVert + 2; 140 ierr = PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);CHKERRQ(ierr); 141 } else { 142 if (numEntries != numFaceVert + 2) { 143 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for mixed faces in Fluent files"); 144 } 145 } 146 } 147 for (e = 0; e < numEntries; e++) { 148 ierr = PetscViewerRead(viewer, buffer, 1, PETSC_STRING);CHKERRQ(ierr); 149 snum = sscanf(buffer, "%x", &(((PetscInt*)s->data)[f*numEntries + e])); 150 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 151 } 152 } 153 s->nd = numEntries - 2; 154 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr); 155 } 156 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr); 157 158 } else { /* Unknown section type */ 159 PetscInt depth = 1; 160 do { 161 /* Match parentheses when parsing unknown sections */ 162 do {ierr = PetscViewerRead(viewer, &(buffer[0]), 1, PETSC_CHAR);CHKERRQ(ierr);} 163 while (buffer[0] != '(' && buffer[0] != ')'); 164 if (buffer[0] == '(') depth++; 165 if (buffer[0] == ')') depth--; 166 } while (depth > 0); 167 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '\n');CHKERRQ(ierr); 168 } 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNCT__ 173 #define __FUNCT__ "DMPlexCreateFluent" 174 /*@C 175 DMPlexCreateFluent - Create a DMPlex mesh from a Fluent mesh file. 176 177 Collective on comm 178 179 Input Parameters: 180 + comm - The MPI communicator 181 . viewer - The Viewer associated with a Fluent mesh file 182 - interpolate - Create faces and edges in the mesh 183 184 Output Parameter: 185 . dm - The DM object representing the mesh 186 187 Note: http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm 188 189 Level: beginner 190 191 .keywords: mesh, fluent, case 192 .seealso: DMPLEX, DMCreate() 193 @*/ 194 PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm) 195 { 196 PetscMPIInt rank; 197 PetscInt c, f, v, dim = PETSC_DETERMINE, numCells = 0, numVertices = 0, numCellVertices = PETSC_DETERMINE; 198 PetscInt numFaces = PETSC_DETERMINE, numFaceEntries = PETSC_DETERMINE, numFaceVertices = PETSC_DETERMINE; 199 PetscInt *faces = NULL, *cellVertices; 200 PetscInt d, coordSize; 201 PetscScalar *coords, *coordsIn = NULL; 202 PetscSection coordSection; 203 Vec coordinates; 204 PetscErrorCode ierr; 205 206 PetscFunctionBegin; 207 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 208 209 if (!rank) { 210 FluentSection s; 211 do { 212 ierr = DMPlexCreateFluent_ReadSection(viewer, &s);CHKERRQ(ierr); 213 if (s.index == 2) { /* Dimension */ 214 dim = s.nd; 215 216 } else if (s.index == 10) { /* Vertices */ 217 if (s.zoneID == 0) numVertices = s.last; 218 else { 219 if (coordsIn) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files"); 220 coordsIn = s.data; 221 } 222 223 } else if (s.index == 12) { /* Cells */ 224 if (s.zoneID == 0) numCells = s.last; 225 else { 226 switch (s.nd) { 227 case 0: numCellVertices = PETSC_DETERMINE; break; 228 case 1: numCellVertices = 3; break; /* triangular */ 229 case 2: numCellVertices = 4; break; /* tetrahedral */ 230 case 3: numCellVertices = 4; break; /* quadrilateral */ 231 case 4: numCellVertices = 8; break; /* hexahedral */ 232 case 5: numCellVertices = 5; break; /* pyramid */ 233 case 6: numCellVertices = 6; break; /* wedge */ 234 default: numCellVertices = PETSC_DETERMINE; 235 } 236 } 237 238 } else if (s.index == 13) { /* Facets */ 239 if (s.zoneID == 0) { /* Header section */ 240 numFaces = s.last - s.first + 1; 241 if (s.nd == 0) numFaceVertices = PETSC_DETERMINE; 242 else numFaceVertices = s.nd; 243 } else { /* Data section */ 244 if (numFaceVertices != PETSC_DETERMINE && s.nd != numFaceVertices) { 245 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed facets in Fluent files are not supported"); 246 } 247 if (numFaces < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file"); 248 if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd; 249 numFaceEntries = numFaceVertices + 2; 250 if (!faces) {ierr = PetscMalloc1(numFaces*numFaceEntries, &faces);CHKERRQ(ierr);} 251 ierr = PetscMemcpy(&(faces[(s.first-1)*numFaceEntries]), s.data, (s.last-s.first+1)*numFaceEntries*sizeof(PetscInt));CHKERRQ(ierr); 252 ierr = PetscFree(s.data);CHKERRQ(ierr); 253 } 254 } 255 } while (s.index >= 0); 256 } 257 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 258 if (dim < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension"); 259 260 /* Allocate cell-vertex mesh */ 261 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 262 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 263 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 264 ierr = DMPlexSetChart(*dm, 0, numCells + numVertices);CHKERRQ(ierr); 265 if (!rank) { 266 if (numCells < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown number of cells in Fluent file"); 267 /* If no cell type was given we assume simplices */ 268 if (numCellVertices == PETSC_DETERMINE) numCellVertices = numFaceVertices + 1; 269 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(*dm, c, numCellVertices);CHKERRQ(ierr);} 270 } 271 ierr = DMSetUp(*dm);CHKERRQ(ierr); 272 273 if (!rank) { 274 /* Derive cell-vertex list from face-vertex and face-cell maps */ 275 ierr = PetscMalloc1(numCells*numCellVertices, &cellVertices);CHKERRQ(ierr); 276 for (c = 0; c < numCells*numCellVertices; c++) cellVertices[c] = -1; 277 for (f = 0; f < numFaces; f++) { 278 PetscInt *cell; 279 const PetscInt cl = faces[f*numFaceEntries + numFaceVertices]; 280 const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1]; 281 const PetscInt *face = &(faces[f*numFaceEntries]); 282 283 if (cl > 0) { 284 cell = &(cellVertices[(cl-1) * numCellVertices]); 285 for (v = 0; v < numFaceVertices; v++) { 286 PetscBool found = PETSC_FALSE; 287 for (c = 0; c < numCellVertices; c++) { 288 if (cell[c] < 0) break; 289 if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;} 290 } 291 if (!found) cell[c] = face[v]-1 + numCells; 292 } 293 } 294 if (cr > 0) { 295 cell = &(cellVertices[(cr-1) * numCellVertices]); 296 for (v = 0; v < numFaceVertices; v++) { 297 PetscBool found = PETSC_FALSE; 298 for (c = 0; c < numCellVertices; c++) { 299 if (cell[c] < 0) break; 300 if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;} 301 } 302 if (!found) cell[c] = face[v]-1 + numCells; 303 } 304 } 305 } 306 for (c = 0; c < numCells; c++) { 307 ierr = DMPlexSetCone(*dm, c, &(cellVertices[c*numCellVertices]));CHKERRQ(ierr); 308 } 309 } 310 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 311 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 312 if (interpolate) { 313 DM idm = NULL; 314 315 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 316 ierr = DMDestroy(dm);CHKERRQ(ierr); 317 *dm = idm; 318 } 319 320 /* Read coordinates */ 321 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 322 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 323 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 324 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 325 for (v = numCells; v < numCells+numVertices; ++v) { 326 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 327 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 328 } 329 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 330 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 331 ierr = VecCreate(comm, &coordinates);CHKERRQ(ierr); 332 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 333 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 334 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 335 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 336 if (!rank) { 337 for (v = 0; v < numVertices; ++v) { 338 for (d = 0; d < dim; ++d) { 339 coords[v*dim+d] = coordsIn[v*dim+d]; 340 } 341 } 342 } 343 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 344 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 345 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 346 if (!rank) { 347 ierr = PetscFree(cellVertices);CHKERRQ(ierr); 348 ierr = PetscFree(faces);CHKERRQ(ierr); 349 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 350 } 351 PetscFunctionReturn(0); 352 } 353