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 } 96 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr); 97 98 } else if (s->index == 13) { /* Faces */ 99 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr); 100 snum = sscanf(buffer, "(%x", &(s->zoneID)); 101 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 102 if (s->zoneID == 0) { /* Header section */ 103 snum = sscanf(buffer, "(%x %x %x %d)", &(s->zoneID), &(s->first), &(s->last), &(s->nd)); 104 if (snum != 4) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 105 } else { /* Data section */ 106 PetscInt f, e, numEntries, numFaces; 107 snum = sscanf(buffer, "(%x %x %x %d %d)", &(s->zoneID), &(s->first), &(s->last), &(s->type), &(s->nd)); 108 if (snum != 5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 109 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '(');CHKERRQ(ierr); 110 switch (s->nd) { 111 case 0: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed faces in Fluent files are not supported"); 112 case 2: numEntries = 2 + 2; break; /* linear */ 113 case 3: numEntries = 2 + 3; break; /* triangular */ 114 case 4: numEntries = 2 + 4; break; /* quadrilateral */ 115 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown face type in Fluent file"); 116 } 117 numFaces = s->last-s->first + 1; 118 ierr = PetscMalloc1(numEntries*numFaces, (PetscInt**)&s->data);CHKERRQ(ierr); 119 for (f = 0; f < numFaces; f++) { 120 for (e = 0; e < numEntries; e++) { 121 ierr = PetscViewerRead(viewer, buffer, 1, PETSC_STRING);CHKERRQ(ierr); 122 snum = sscanf(buffer, "%x", &(((PetscInt*)s->data)[f*numEntries + e])); 123 if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "File is not a valid Fluent file"); 124 } 125 } 126 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr); 127 } 128 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, ')');CHKERRQ(ierr); 129 130 } else { /* Unknown section type */ 131 PetscInt depth = 1; 132 do { 133 /* Match parentheses when parsing unknown sections */ 134 do {ierr = PetscViewerRead(viewer, &(buffer[0]), 1, PETSC_CHAR);CHKERRQ(ierr);} 135 while (buffer[0] != '(' && buffer[0] != ')'); 136 if (buffer[0] == '(') depth++; 137 if (buffer[0] == ')') depth--; 138 } while (depth > 0); 139 ierr = DMPlexCreateFluent_ReadString(viewer, buffer, '\n');CHKERRQ(ierr); 140 } 141 PetscFunctionReturn(0); 142 } 143 144 #undef __FUNCT__ 145 #define __FUNCT__ "DMPlexCreateFluent" 146 /*@C 147 DMPlexCreateFluent - Create a DMPlex mesh from a Fluent mesh file. 148 149 Collective on comm 150 151 Input Parameters: 152 + comm - The MPI communicator 153 . viewer - The Viewer associated with a Fluent mesh file 154 - interpolate - Create faces and edges in the mesh 155 156 Output Parameter: 157 . dm - The DM object representing the mesh 158 159 Note: http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1490.htm 160 161 Level: beginner 162 163 .keywords: mesh, fluent, case 164 .seealso: DMPLEX, DMCreate() 165 @*/ 166 PetscErrorCode DMPlexCreateFluent(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm) 167 { 168 PetscMPIInt rank; 169 PetscInt c, f, v, dim = PETSC_DETERMINE, numCells = 0, numVertices = 0, numCellVertices = PETSC_DETERMINE; 170 PetscInt numFaces = PETSC_DETERMINE, numFaceEntries = PETSC_DETERMINE, numFaceVertices = PETSC_DETERMINE; 171 PetscInt *faces = NULL, *cellVertices; 172 PetscInt d, coordSize; 173 PetscScalar *coords, *coordsIn = NULL; 174 PetscSection coordSection; 175 Vec coordinates; 176 PetscErrorCode ierr; 177 178 PetscFunctionBegin; 179 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 180 181 if (!rank) { 182 FluentSection s; 183 do { 184 ierr = DMPlexCreateFluent_ReadSection(viewer, &s);CHKERRQ(ierr); 185 if (s.index == 2) { /* Dimension */ 186 dim = s.nd; 187 188 } else if (s.index == 10) { /* Vertices */ 189 if (s.zoneID == 0) numVertices = s.last; 190 else { 191 if (coordsIn) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Currently no support for multiple coordinate sets in Fluent files"); 192 coordsIn = s.data; 193 } 194 195 } else if (s.index == 12) { /* Cells */ 196 if (s.zoneID == 0) numCells = s.last; 197 else { 198 switch (s.nd) { 199 case 0: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed elements in Fluent files are not supported"); 200 case 1: numCellVertices = 3; break; /* triangular */ 201 case 2: numCellVertices = 4; break; /* tetrahedral */ 202 case 3: numCellVertices = 4; break; /* quadrilateral */ 203 case 4: numCellVertices = 8; break; /* hexahedral */ 204 case 5: numCellVertices = 5; break; /* pyramid */ 205 case 6: numCellVertices = 6; break; /* wedge */ 206 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown cell element-type in Fluent file"); 207 } 208 } 209 210 } else if (s.index == 13) { /* Facets */ 211 if (s.zoneID == 0) { /* Header section */ 212 numFaces = s.last - s.first + 1; 213 } else { /* Data section */ 214 if (s.nd == 0 || (numFaceVertices != PETSC_DETERMINE && s.nd != numFaceVertices)) { 215 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mixed facets in Fluent files are currently not supported"); 216 } 217 if (numFaces == PETSC_DETERMINE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No header section for facets in Fluent file"); 218 if (numFaceVertices == PETSC_DETERMINE) numFaceVertices = s.nd; 219 numFaceEntries = numFaceVertices + 2; 220 if (!faces) {ierr = PetscMalloc1(numFaces*numFaceEntries, &faces);CHKERRQ(ierr);} 221 ierr = PetscMemcpy(&(faces[(s.first-1)*numFaceEntries]), s.data, (s.last-s.first+1)*numFaceEntries*sizeof(PetscInt));CHKERRQ(ierr); 222 ierr = PetscFree(s.data);CHKERRQ(ierr); 223 } 224 } 225 } while (s.index >= 0); 226 } 227 ierr = MPI_Bcast(&dim, 1, MPIU_INT, 0, comm);CHKERRQ(ierr); 228 if (dim < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Fluent file does not include dimension"); 229 230 /* Allocate cell-vertex mesh */ 231 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 232 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 233 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr); 234 ierr = DMPlexSetChart(*dm, 0, numCells + numVertices);CHKERRQ(ierr); 235 if (!rank) { 236 if (numCellVertices == PETSC_DETERMINE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cell type not specified in Fluent file"); 237 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(*dm, c, numCellVertices);CHKERRQ(ierr);} 238 } 239 ierr = DMSetUp(*dm);CHKERRQ(ierr); 240 241 if (!rank) { 242 /* Derive cell-vertex list from face-vertex and face-cell maps */ 243 if (numCells == PETSC_DETERMINE || numCellVertices == PETSC_DETERMINE) { 244 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Insufficent cell header information in Fluent file"); 245 } 246 ierr = PetscMalloc1(numCells*numCellVertices, &cellVertices);CHKERRQ(ierr); 247 for (c = 0; c < numCells*numCellVertices; c++) cellVertices[c] = -1; 248 for (f = 0; f < numFaces; f++) { 249 PetscInt *cell; 250 const PetscInt cl = faces[f*numFaceEntries + numFaceVertices]; 251 const PetscInt cr = faces[f*numFaceEntries + numFaceVertices + 1]; 252 const PetscInt *face = &(faces[f*numFaceEntries]); 253 254 if (cl > 0) { 255 cell = &(cellVertices[(cl-1) * numCellVertices]); 256 for (v = 0; v < numFaceVertices; v++) { 257 PetscBool found = PETSC_FALSE; 258 for (c = 0; c < numCellVertices; c++) { 259 if (cell[c] < 0) break; 260 if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;} 261 } 262 if (!found) cell[c] = face[v]-1 + numCells; 263 } 264 } 265 if (cr > 0) { 266 cell = &(cellVertices[(cr-1) * numCellVertices]); 267 for (v = 0; v < numFaceVertices; v++) { 268 PetscBool found = PETSC_FALSE; 269 for (c = 0; c < numCellVertices; c++) { 270 if (cell[c] < 0) break; 271 if (cell[c] == face[v]-1 + numCells) {found = PETSC_TRUE; break;} 272 } 273 if (!found) cell[c] = face[v]-1 + numCells; 274 } 275 } 276 } 277 for (c = 0; c < numCells; c++) { 278 ierr = DMPlexSetCone(*dm, c, &(cellVertices[c*numCellVertices]));CHKERRQ(ierr); 279 } 280 } 281 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr); 282 ierr = DMPlexStratify(*dm);CHKERRQ(ierr); 283 if (interpolate) { 284 DM idm = NULL; 285 286 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr); 287 ierr = DMDestroy(dm);CHKERRQ(ierr); 288 *dm = idm; 289 } 290 291 /* Read coordinates */ 292 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 293 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr); 294 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr); 295 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr); 296 for (v = numCells; v < numCells+numVertices; ++v) { 297 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr); 298 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr); 299 } 300 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr); 301 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr); 302 ierr = VecCreate(comm, &coordinates);CHKERRQ(ierr); 303 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr); 304 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr); 305 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr); 306 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 307 if (!rank) { 308 for (v = 0; v < numVertices; ++v) { 309 for (d = 0; d < dim; ++d) { 310 coords[v*dim+d] = coordsIn[v*dim+d]; 311 } 312 } 313 } 314 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 315 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr); 316 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 317 if (!rank) { 318 ierr = PetscFree(cellVertices);CHKERRQ(ierr); 319 ierr = PetscFree(faces);CHKERRQ(ierr); 320 ierr = PetscFree(coordsIn);CHKERRQ(ierr); 321 } 322 PetscFunctionReturn(0); 323 } 324