FreeFOAM The Cross-Platform CFD Toolkit
netgenNeutralToFoam.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Application
25  netgenNeutralToFoam
26 
27 Description
28  Converts neutral file format as written by Netgen v4.4.
29 
30  Example:
31 
32  9
33  1.000000 1.000000 1.000000
34  0.000000 1.000000 1.000000
35  0.000000 0.000000 1.000000
36  1.000000 0.000000 1.000000
37  0.000000 1.000000 0.000000
38  1.000000 1.000000 0.000000
39  1.000000 0.000000 0.000000
40  0.000000 0.000000 0.000000
41  0.500000 0.500000 0.500000
42  12
43  1 7 8 9 3
44  1 5 9 6 8
45  1 5 9 2 1
46  1 4 9 7 6
47  1 7 8 6 9
48  1 4 6 1 9
49  1 5 9 8 2
50  1 4 1 2 9
51  1 1 6 5 9
52  1 2 3 4 9
53  1 8 9 3 2
54  1 4 9 3 7
55  12
56  1 1 2 4
57  1 3 4 2
58  2 5 6 8
59  2 7 8 6
60  3 1 4 6
61  3 7 6 4
62  5 2 1 5
63  5 6 5 1
64  5 3 2 8
65  5 5 8 2
66  6 4 3 7
67  6 8 7 3
68 
69 NOTE:
70  - reverse order of boundary faces using geometric test.
71  (not very space efficient)
72  - order of tet vertices only tested on one file.
73  - all patch/cell/vertex numbers offset by one.
74 
75 Usage
76 
77  - netgenNeutralToFoam [OPTIONS] <Neutral file>
78 
79  @param <Neutral file> \n
80  @todo Detailed description of argument.
81 
82  @param -case <dir>\n
83  Case directory.
84 
85  @param -parallel \n
86  Run in parallel.
87 
88  @param -help \n
89  Display help message.
90 
91  @param -doc \n
92  Display Doxygen API documentation page for this application.
93 
94  @param -srcDoc \n
95  Display Doxygen source documentation page for this application.
96 
97 \*---------------------------------------------------------------------------*/
98 
99 #include <OpenFOAM/argList.H>
100 #include <OpenFOAM/Time.H>
101 #include <OpenFOAM/polyMesh.H>
102 #include <OpenFOAM/IFstream.H>
103 #include <OpenFOAM/polyPatch.H>
104 #include <OpenFOAM/cellModeller.H>
105 #include <OpenFOAM/triFace.H>
106 
107 using namespace Foam;
108 
109 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
110 
111 // Main program:
112 
113 int main(int argc, char *argv[])
114 {
115  argList::validArgs.append("Neutral file");
116 
117 # include <OpenFOAM/setRootCase.H>
118 # include <OpenFOAM/createTime.H>
119 
120  fileName neuFile(args.additionalArgs()[0]);
121 
122 
123  IFstream str(neuFile);
124 
125 
126  //
127  // Read nodes.
128  //
129  label nNodes(readLabel(str));
130 
131  Info<< "nNodes:" << nNodes << endl;
132 
133 
134  pointField points(nNodes);
135 
136  forAll(points, pointI)
137  {
138  scalar x,y,z;
139 
140  str >> x >> y >> z;
141 
142  points[pointI] = point(x, y, z);
143  }
144 
145 
146 
147 
148  label nTets(readLabel(str));
149 
150  Info<< "nTets:" << nTets << endl;
151 
152  const cellModel& tet = *(cellModeller::lookup("tet"));
153 
154  cellShapeList cells(nTets);
155 
156  labelList tetPoints(4);
157 
158  forAll(cells, cellI)
159  {
160  label domain(readLabel(str));
161 
162  if (domain != 1)
163  {
165  << "Cannot handle multiple domains"
166  << nl << "Ignoring domain " << domain << " setting on line "
167  << str.lineNumber() << endl;
168  }
169 
170  tetPoints[1] = readLabel(str) - 1;
171  tetPoints[0] = readLabel(str) - 1;
172  tetPoints[2] = readLabel(str) - 1;
173  tetPoints[3] = readLabel(str) - 1;
174 
175  cells[cellI] = cellShape(tet, tetPoints);
176  }
177 
178 
179 
180  label nFaces(readLabel(str));
181 
182  Info<< "nFaces:" << nFaces << endl;
183 
184  // Unsorted boundary faces
185  faceList boundaryFaces(nFaces);
186 
187  // For each boundary faces the Foam patchID
188  labelList boundaryPatch(nFaces, -1);
189 
190  // Max patch.
191  label maxPatch = 0;
192 
193  // Boundary faces as three vertices
194  HashTable<label, triFace, Hash<triFace> > vertsToBoundary(nFaces);
195 
196  forAll(boundaryFaces, faceI)
197  {
198  label patchI(readLabel(str));
199 
200  if (patchI < 0)
201  {
203  << "Invalid boundary region number " << patchI
204  << " on line " << str.lineNumber()
205  << exit(FatalError);
206  }
207 
208 
209  maxPatch = max(maxPatch, patchI);
210 
211  triFace tri(readLabel(str)-1, readLabel(str)-1, readLabel(str)-1);
212 
213  // Store boundary face as is for now. Later on reverse it.
214  boundaryFaces[faceI].setSize(3);
215  boundaryFaces[faceI][0] = tri[0];
216  boundaryFaces[faceI][1] = tri[1];
217  boundaryFaces[faceI][2] = tri[2];
218  boundaryPatch[faceI] = patchI;
219 
220  vertsToBoundary.insert(tri, faceI);
221  }
222 
223  label nPatches = maxPatch + 1;
224 
225 
226  // Use hash of points to get owner cell and orient the boundary face.
227  // For storage reasons I store the triangles and loop over the cells instead
228  // of the other way around (store cells and loop over triangles) though
229  // that would be faster.
230  forAll(cells, cellI)
231  {
232  const cellShape& cll = cells[cellI];
233 
234  // Get the four (outwards pointing) faces of the cell
235  faceList tris(cll.faces());
236 
237  forAll(tris, i)
238  {
239  const face& f = tris[i];
240 
241  // Is there any boundary face with same vertices?
242  // (uses commutative hash)
243  HashTable<label, triFace, Hash<triFace> >::iterator iter =
244  vertsToBoundary.find(triFace(f[0], f[1], f[2]));
245 
246  if (iter != vertsToBoundary.end())
247  {
248  label faceI = iter();
249  const triFace& tri = iter.key();
250 
251  // Determine orientation of tri v.s. cell centre.
252  point cc(cll.centre(points));
253  point fc(tri.centre(points));
254  vector fn(tri.normal(points));
255 
256  if (((fc - cc) & fn) < 0)
257  {
258  // Boundary face points inwards. Flip.
259  boundaryFaces[faceI] = boundaryFaces[faceI].reverseFace();
260  }
261 
262  // Done this face so erase from hash
263  vertsToBoundary.erase(iter);
264  }
265  }
266  }
267 
268 
269  if (vertsToBoundary.size())
270  {
271  // Didn't find cells connected to boundary faces.
273  << "There are boundary faces without attached cells."
274  << "Boundary faces (as triFaces):" << vertsToBoundary.toc()
275  << endl;
276  }
277 
278 
279  // Storage for boundary faces sorted into patches
280 
281  faceListList patchFaces(nPatches);
282 
283  wordList patchNames(nPatches);
284 
285  forAll(patchNames, patchI)
286  {
287  patchNames[patchI] = word("patch") + name(patchI);
288  }
289 
290  wordList patchTypes(nPatches, polyPatch::typeName);
291  word defaultFacesName = "defaultFaces";
292  word defaultFacesType = polyPatch::typeName;
293  wordList patchPhysicalTypes(nPatches, polyPatch::typeName);
294 
295  {
296  // Sort boundaryFaces by patch.
297  List<DynamicList<face> > allPatchFaces(nPatches);
298 
299  forAll(boundaryPatch, faceI)
300  {
301  label patchI = boundaryPatch[faceI];
302 
303  allPatchFaces[patchI].append(boundaryFaces[faceI]);
304  }
305 
306  Info<< "Patches:" << nl
307  << "\tNeutral Boundary\tPatch name\tSize" << nl
308  << "\t----------------\t----------\t----" << endl;
309 
310  forAll(allPatchFaces, patchI)
311  {
312  Info<< '\t' << patchI << "\t\t\t"
313  << patchNames[patchI] << "\t\t"
314  << allPatchFaces[patchI].size() << endl;
315 
316  patchFaces[patchI].transfer(allPatchFaces[patchI]);
317  }
318 
319  Info<< endl;
320  }
321 
322 
323  polyMesh mesh
324  (
325  IOobject
326  (
328  runTime.constant(),
329  runTime
330  ),
331  xferMove(points),
332  cells,
333  patchFaces,
334  patchNames,
335  patchTypes,
339  );
340 
341  Info<< "Writing mesh to " << runTime.constant() << endl << endl;
342 
343  mesh.write();
344 
345 
346  Info<< "End\n" << endl;
347 
348  return 0;
349 }
350 
351 
352 // ************************ vim: set sw=4 sts=4 et: ************************ //