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