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) 2011-2016 OpenFOAM Foundation
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 \*---------------------------------------------------------------------------*/
76 
77 #include "argList.H"
78 #include "Time.H"
79 #include "polyMesh.H"
80 #include "IFstream.H"
81 #include "polyPatch.H"
82 #include "cellModeller.H"
83 #include "triFace.H"
84 
85 using namespace Foam;
86 
87 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
88 
89 
90 int main(int argc, char *argv[])
91 {
92  argList::validArgs.append("Neutral file");
93 
94  #include "setRootCase.H"
95  #include "createTime.H"
96 
97  IFstream str(args[1]);
98 
99  //
100  // Read nodes.
101  //
102  label nNodes(readLabel(str));
103 
104  Info<< "nNodes:" << nNodes << endl;
105 
106  pointField points(nNodes);
107 
108  forAll(points, pointi)
109  {
110  scalar x,y,z;
111 
112  str >> x >> y >> z;
113 
114  points[pointi] = point(x, y, z);
115  }
116 
117 
118 
119 
120  label nTets(readLabel(str));
121 
122  Info<< "nTets:" << nTets << endl;
123 
124  const cellModel& tet = *(cellModeller::lookup("tet"));
125 
126  cellShapeList cells(nTets);
127 
128  labelList tetPoints(4);
129 
130  forAll(cells, celli)
131  {
132  label domain(readLabel(str));
133 
134  if (domain != 1)
135  {
137  << "Cannot handle multiple domains"
138  << nl << "Ignoring domain " << domain << " setting on line "
139  << str.lineNumber() << endl;
140  }
141 
142  tetPoints[1] = readLabel(str) - 1;
143  tetPoints[0] = readLabel(str) - 1;
144  tetPoints[2] = readLabel(str) - 1;
145  tetPoints[3] = readLabel(str) - 1;
146 
147  cells[celli] = cellShape(tet, tetPoints);
148  }
149 
150 
151 
152  label nFaces(readLabel(str));
153 
154  Info<< "nFaces:" << nFaces << endl;
155 
156  // Unsorted boundary faces
157  faceList boundaryFaces(nFaces);
158 
159  // For each boundary faces the Foam patchID
160  labelList boundaryPatch(nFaces, -1);
161 
162  // Max patch.
163  label maxPatch = 0;
164 
165  // Boundary faces as three vertices
166  HashTable<label, triFace, Hash<triFace>> vertsToBoundary(nFaces);
167 
168  forAll(boundaryFaces, facei)
169  {
170  label patchi(readLabel(str));
171 
172  if (patchi < 0)
173  {
175  << "Invalid boundary region number " << patchi
176  << " on line " << str.lineNumber()
177  << exit(FatalError);
178  }
179 
180 
181  maxPatch = max(maxPatch, patchi);
182 
183  triFace tri(readLabel(str)-1, readLabel(str)-1, readLabel(str)-1);
184 
185  // Store boundary face as is for now. Later on reverse it.
186  boundaryFaces[facei].setSize(3);
187  boundaryFaces[facei][0] = tri[0];
188  boundaryFaces[facei][1] = tri[1];
189  boundaryFaces[facei][2] = tri[2];
190  boundaryPatch[facei] = patchi;
191 
192  vertsToBoundary.insert(tri, facei);
193  }
194 
195  label nPatches = maxPatch + 1;
196 
197 
198  // Use hash of points to get owner cell and orient the boundary face.
199  // For storage reasons I store the triangles and loop over the cells instead
200  // of the other way around (store cells and loop over triangles) though
201  // that would be faster.
202  forAll(cells, celli)
203  {
204  const cellShape& cll = cells[celli];
205 
206  // Get the four (outwards pointing) faces of the cell
207  faceList tris(cll.faces());
208 
209  forAll(tris, i)
210  {
211  const face& f = tris[i];
212 
213  // Is there any boundary face with same vertices?
214  // (uses commutative hash)
216  vertsToBoundary.find(triFace(f[0], f[1], f[2]));
217 
218  if (iter != vertsToBoundary.end())
219  {
220  label facei = iter();
221  const triFace& tri = iter.key();
222 
223  // Determine orientation of tri v.s. cell centre.
224  point cc(cll.centre(points));
225  point fc(tri.centre(points));
226  vector fn(tri.normal(points));
227 
228  if (((fc - cc) & fn) < 0)
229  {
230  // Boundary face points inwards. Flip.
231  boundaryFaces[facei].flip();
232  }
233 
234  // Done this face so erase from hash
235  vertsToBoundary.erase(iter);
236  }
237  }
238  }
239 
240 
241  if (vertsToBoundary.size())
242  {
243  // Didn't find cells connected to boundary faces.
245  << "There are boundary faces without attached cells."
246  << "Boundary faces (as triFaces):" << vertsToBoundary.toc()
247  << endl;
248  }
249 
250 
251  // Storage for boundary faces sorted into patches
252 
253  faceListList patchFaces(nPatches);
254 
255  wordList patchNames(nPatches);
256 
257  forAll(patchNames, patchi)
258  {
259  patchNames[patchi] = word("patch") + name(patchi);
260  }
261 
262  wordList patchTypes(nPatches, polyPatch::typeName);
263  word defaultFacesName = "defaultFaces";
264  word defaultFacesType = polyPatch::typeName;
265  wordList patchPhysicalTypes(nPatches, polyPatch::typeName);
266 
267  {
268  // Sort boundaryFaces by patch.
269  List<DynamicList<face>> allPatchFaces(nPatches);
270 
271  forAll(boundaryPatch, facei)
272  {
273  label patchi = boundaryPatch[facei];
274 
275  allPatchFaces[patchi].append(boundaryFaces[facei]);
276  }
277 
278  Info<< "Patches:" << nl
279  << "\tNeutral Boundary\tPatch name\tSize" << nl
280  << "\t----------------\t----------\t----" << endl;
281 
282  forAll(allPatchFaces, patchi)
283  {
284  Info<< '\t' << patchi << "\t\t\t"
285  << patchNames[patchi] << "\t\t"
286  << allPatchFaces[patchi].size() << endl;
287 
288  patchFaces[patchi].transfer(allPatchFaces[patchi]);
289  }
290 
291  Info<< endl;
292  }
293 
294 
295  polyMesh mesh
296  (
297  IOobject
298  (
300  runTime.constant(),
301  runTime
302  ),
303  xferMove(points),
304  cells,
305  patchFaces,
306  patchNames,
307  patchTypes,
310  patchPhysicalTypes
311  );
312 
313  Info<< "Writing mesh to " << runTime.constant() << endl << endl;
314 
315  mesh.write();
316 
317 
318  Info<< "End\n" << endl;
319 
320  return 0;
321 }
322 
323 
324 // ************************************************************************* //
word defaultFacesType
Definition: readKivaGrid.H:461
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
static const cellModel * lookup(const word &)
Look up a model by name and return a pointer to the model or nullptr.
Definition: cellModeller.C:88
Like polyPatch but without reference to mesh. patchIdentifier::index is not used. Used in boundaryMes...
Definition: boundaryPatch.H:57
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:110
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
virtual bool write(const bool valid=true) const
Write mesh using IO settings from time.
Definition: fvMesh.C:877
An analytical geometric cellShape.
Definition: cellShape.H:69
wordList patchTypes(nPatches)
static word defaultRegion
Return the default region name.
Definition: polyMesh.H:309
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
point centre(const pointField &) const
Return centre (centroid)
Definition: triFaceI.H:163
static SLList< string > validArgs
A list of valid (mandatory) arguments.
Definition: argList.H:154
void transfer(HashTable< T, Key, Hash > &)
Transfer the contents of the argument table into this table.
Definition: HashTable.C:513
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:142
scalar y
faceList faces() const
Faces of this cell.
Definition: cellShapeI.H:180
dynamicFvMesh & mesh
const cellShapeList & cells
const pointField & points
Xfer< T > xferMove(T &)
Construct by transferring the contents of the arg.
A class for handling words, derived from string.
Definition: word.H:59
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:68
face triFace(3)
wordList patchNames(nPatches)
An STL-conforming hash table.
Definition: HashTable.H:62
label readLabel(Istream &is)
Definition: label.H:64
static const char nl
Definition: Ostream.H:262
Input from file stream.
Definition: IFstream.H:81
word defaultFacesName
Definition: readKivaGrid.H:460
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
label patchi
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
vector normal(const pointField &) const
Vector normal; magnitude is equal to area of face.
Definition: triFaceI.H:181
point centre(const pointField &) const
Centroid of the cell.
Definition: cellShapeI.H:259
Maps a geometry to a set of cell primitives, which enables geometric cell data to be calculated witho...
Definition: cellModel.H:64
messageStream Info
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Foam::argList args(argc, argv)
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Namespace for OpenFOAM.