netgenNeutralToFoam.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2019 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.area(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 
256 
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  move(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 // ************************************************************************* //
scalar y
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
An STL-conforming hash table.
Definition: HashTable.H:127
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:167
Input from file stream.
Definition: IFstream.H:85
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
virtual Ostream & write(const char)=0
Write character.
static SLList< string > validArgs
A list of valid (mandatory) arguments.
Definition: argList.H:153
Maps a geometry to a set of cell primitives, which enables geometric cell data to be calculated witho...
Definition: cellModel.H:65
static const cellModel * lookup(const word &)
Look up a model by name and return a pointer to the model or nullptr.
Definition: cellModeller.C:100
An analytical geometric cellShape.
Definition: cellShape.H:72
point centre(const pointField &) const
Centroid of the cell.
Definition: cellShapeI.H:259
faceList faces() const
Faces of this cell.
Definition: cellShapeI.H:180
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:76
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
static word defaultRegion
Return the default region name.
Definition: polyMesh.H:267
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:71
point centre(const pointField &) const
Return centre (centroid)
Definition: triFaceI.H:163
vector area(const pointField &) const
Return vector area.
Definition: triFaceI.H:174
A class for handling words, derived from string.
Definition: word.H:62
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
int main(int argc, char *argv[])
Definition: financialFoam.C:44
label patchi
const pointField & points
const cellShapeList & cells
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
messageStream Info
vector point
Point is a vector.
Definition: point.H:41
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
label readLabel(Istream &is)
Definition: label.H:64
error FatalError
static const char nl
Definition: Ostream.H:266
wordList patchTypes(nPatches)
wordList patchNames(nPatches)
face triFace(3)
labelList f(nPoints)
word defaultFacesName
Definition: readKivaGrid.H:454
word defaultFacesType
Definition: readKivaGrid.H:455
label nPatches
Definition: readKivaGrid.H:396
Foam::argList args(argc, argv)