FreeFOAM The Cross-Platform CFD Toolkit
tetgenToFoam.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  tetgenToFoam
26 
27 Description
28  Converts .ele and .node and .face files, written by tetgen.
29 
30  Make sure to use add boundary attributes to the smesh file
31  (5 fifth column in the element section)
32  and run tetgen with -f option.
33 
34  Sample smesh file:
35  @verbatim
36  # cube.smesh -- A 10x10x10 cube
37  8 3
38  1 0 0 0
39  2 0 10 0
40  3 10 10 0
41  4 10 0 0
42  5 0 0 10
43  6 0 10 10
44  7 10 10 10
45  8 10 0 10
46  6 1 # 1 for boundary info present
47  4 1 2 3 4 11 # region number 11
48  4 5 6 7 8 21 # region number 21
49  4 1 2 6 5 3
50  4 4 3 7 8 43
51  4 1 5 8 4 5
52  4 2 6 7 3 65
53  0
54  0
55  @endverbatim
56 
57 NOTE:
58  - for some reason boundary faces point inwards. I just reverse them
59  always. Might use some geometric check instead.
60  - marked faces might not actually be boundary faces of mesh.
61  This is hopefully handled now by first constructing without boundaries
62  and then reconstructing with boundary faces.
63 
64 Usage
65 
66  - tetgenToFoam [OPTIONS] <file prefix>
67 
68  @param <file prefix> \n
69  @todo Detailed description of argument.
70 
71  @param -noFaceFile \n
72  Ignore the face file.
73 
74  @param -case <dir>\n
75  Case directory.
76 
77  @param -parallel \n
78  Run in parallel.
79 
80  @param -help \n
81  Display help message.
82 
83  @param -doc \n
84  Display Doxygen API documentation page for this application.
85 
86  @param -srcDoc \n
87  Display Doxygen source documentation page for this application.
88 
89 \*---------------------------------------------------------------------------*/
90 
91 #include <OpenFOAM/argList.H>
92 #include <OpenFOAM/Time.H>
93 #include <OpenFOAM/polyMesh.H>
94 #include <OpenFOAM/IFstream.H>
95 #include <OpenFOAM/cellModeller.H>
96 
97 using namespace Foam;
98 
99 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
100 
101 // Find label of face.
102 label findFace(const primitiveMesh& mesh, const face& f)
103 {
104  const labelList& pFaces = mesh.pointFaces()[f[0]];
105 
106  forAll(pFaces, i)
107  {
108  label faceI = pFaces[i];
109 
110  if (mesh.faces()[faceI] == f)
111  {
112  return faceI;
113  }
114  }
115 
116  FatalErrorIn("findFace(const primitiveMesh&, const face&)")
117  << "Cannot find face " << f << " in mesh." << abort(FatalError);
118 
119  return -1;
120 }
121 
122 
123 // Main program:
124 
125 int main(int argc, char *argv[])
126 {
127  argList::validArgs.append("file prefix");
128  argList::validOptions.insert("noFaceFile", "");
129 
130 # include <OpenFOAM/setRootCase.H>
131 # include <OpenFOAM/createTime.H>
132 
133 
134  bool readFaceFile = !args.optionFound("noFaceFile");
135 
136  fileName prefix(args.additionalArgs()[0]);
137 
138  fileName nodeFile(prefix + ".node");
139  fileName eleFile(prefix + ".ele");
140  fileName faceFile(prefix + ".face");
141 
142  if (!readFaceFile)
143  {
144  Info<< "Files:" << endl
145  << " nodes : " << nodeFile << endl
146  << " elems : " << eleFile << endl
147  << endl;
148  }
149  else
150  {
151  Info<< "Files:" << endl
152  << " nodes : " << nodeFile << endl
153  << " elems : " << eleFile << endl
154  << " faces : " << faceFile << endl
155  << endl;
156 
157  Info<< "Reading .face file for boundary information" << nl << endl;
158  }
159 
160  if (!isFile(nodeFile) || !isFile(eleFile))
161  {
163  << "Cannot read " << nodeFile << " or " << eleFile
164  << exit(FatalError);
165  }
166 
167  if (readFaceFile && !isFile(faceFile))
168  {
170  << "Cannot read " << faceFile << endl
171  << "Did you run tetgen with -f option?" << endl
172  << "If you don't want to read the .face file and thus not have"
173  << " patches please\nrerun with the -noFaceFile option"
174  << exit(FatalError);
175  }
176 
177 
178  IFstream nodeStream(nodeFile);
179 
180  //
181  // Read nodes.
182  //
183 
184  // Read header.
185  string line;
186 
187  do
188  {
189  nodeStream.getLine(line);
190  }
191  while (line.size() && line[0] == '#');
192 
193  IStringStream nodeLine(line);
194 
195  label nNodes, nDims, nNodeAttr;
196  bool hasRegion;
197 
198  nodeLine >> nNodes >> nDims >> nNodeAttr >> hasRegion;
199 
200 
201  Info<< "Read .node header:" << endl
202  << " nodes : " << nNodes << endl
203  << " nDims : " << nDims << endl
204  << " nAttr : " << nNodeAttr << endl
205  << " hasRegion : " << hasRegion << endl
206  << endl;
207 
208  //
209  // read points
210  //
211 
212  pointField points(nNodes);
213  Map<label> nodeToPoint(nNodes);
214 
215  {
216  labelList pointIndex(nNodes);
217 
218  label pointI = 0;
219 
220  while (nodeStream.good())
221  {
222  nodeStream.getLine(line);
223 
224  if (line.size() && line[0] != '#')
225  {
226  IStringStream nodeLine(line);
227 
228  label nodeI;
229  scalar x, y, z;
230  label dummy;
231 
232  nodeLine >> nodeI >> x >> y >> z;
233 
234  for (label i = 0; i < nNodeAttr; i++)
235  {
236  nodeLine >> dummy;
237  }
238 
239  if (hasRegion)
240  {
241  nodeLine >> dummy;
242  }
243 
244  // Store point and node number.
245  points[pointI] = point(x, y, z);
246  nodeToPoint.insert(nodeI, pointI);
247  pointI++;
248  }
249  }
250  if (pointI != nNodes)
251  {
252  FatalIOErrorIn(args.executable().c_str(), nodeStream)
253  << "Only " << pointI << " nodes present instead of " << nNodes
254  << " from header." << exit(FatalIOError);
255  }
256  }
257 
258  //
259  // read elements
260  //
261 
262  IFstream eleStream(eleFile);
263 
264  do
265  {
266  eleStream.getLine(line);
267  }
268  while (line.size() && line[0] == '#');
269 
270  IStringStream eleLine(line);
271 
272  label nTets, nPtsPerTet, nElemAttr;
273 
274  eleLine >> nTets >> nPtsPerTet >> nElemAttr;
275 
276 
277  Info<< "Read .ele header:" << endl
278  << " tets : " << nTets << endl
279  << " pointsPerTet : " << nPtsPerTet << endl
280  << " nAttr : " << nElemAttr << endl
281  << endl;
282 
283  if (nPtsPerTet != 4)
284  {
285  FatalIOErrorIn(args.executable().c_str(), eleStream)
286  << "Cannot handle tets with "
287  << nPtsPerTet << " points per tetrahedron in .ele file" << endl
288  << "Can only handle tetrahedra with four points"
289  << exit(FatalIOError);
290  }
291 
292  if (nElemAttr != 0)
293  {
295  << "Element attributes (third elemenent in .ele header)"
296  << " not used" << endl;
297  }
298 
299 
300 
301  const cellModel& tet = *(cellModeller::lookup("tet"));
302 
303  labelList tetPoints(4);
304 
305  cellShapeList cells(nTets);
306  label cellI = 0;
307 
308  while (eleStream.good())
309  {
310  eleStream.getLine(line);
311 
312  if (line.size() && line[0] != '#')
313  {
314  IStringStream eleLine(line);
315 
316  label elemI;
317  eleLine >> elemI;
318 
319  for (label i = 0; i < 4; i++)
320  {
321  label nodeI;
322  eleLine >> nodeI;
323  tetPoints[i] = nodeToPoint[nodeI];
324  }
325 
326  cells[cellI++] = cellShape(tet, tetPoints);
327 
328  // Skip attributes
329  for (label i = 0; i < nElemAttr; i++)
330  {
331  label dummy;
332 
333  eleLine >> dummy;
334  }
335  }
336  }
337 
338 
339  //
340  // Construct mesh with default boundary only
341  //
342 
344  (
345  new polyMesh
346  (
347  IOobject
348  (
350  runTime.constant(),
351  runTime
352  ),
353  xferCopy(points),
354  cells,
355  faceListList(0),
356  wordList(0), // boundaryPatchNames
357  wordList(0), // boundaryPatchTypes
358  "defaultFaces",
359  polyPatch::typeName,
360  wordList(0)
361  )
362  );
363  const polyMesh& mesh = meshPtr;
364 
365 
366 
367  if (readFaceFile)
368  {
369  label nPatches = 0;
370 
371  // List of Foam vertices per boundary face
372  faceList boundaryFaces;
373 
374  // For each boundary faces the Foam patchID
376 
377  //
378  // read boundary faces
379  //
380 
381  IFstream faceStream(faceFile);
382 
383  do
384  {
385  faceStream.getLine(line);
386  }
387  while (line.size() && line[0] == '#');
388 
389  IStringStream faceLine(line);
390 
391  label nFaces, nFaceAttr;
392 
393  faceLine >> nFaces >> nFaceAttr;
394 
395 
396  Info<< "Read .face header:" << endl
397  << " faces : " << nFaces << endl
398  << " nAttr : " << nFaceAttr << endl
399  << endl;
400 
401 
402  if (nFaceAttr != 1)
403  {
404  FatalIOErrorIn(args.executable().c_str(), faceStream)
405  << "Expect boundary markers to be"
406  << " present in .face file." << endl
407  << "This is the second number in the header which is now:"
408  << nFaceAttr << exit(FatalIOError);
409  }
410 
411  // List of Foam vertices per boundary face
412  boundaryFaces.setSize(nFaces);
413 
414  // For each boundary faces the Foam patchID
415  boundaryPatch.setSize(nFaces);
416  boundaryPatch = -1;
417 
418  label faceI = 0;
419 
420  // Region to patch conversion
421  Map<label> regionToPatch;
422 
423  face f(3);
424 
425  while (faceStream.good())
426  {
427  faceStream.getLine(line);
428 
429  if (line.size() && line[0] != '#')
430  {
431  IStringStream faceLine(line);
432 
433  label tetGenFaceI, dummy, region;
434 
435  faceLine >> tetGenFaceI;
436 
437  // Read face and reverse orientation (Foam needs outwards
438  // pointing)
439  for (label i = 0; i < 3; i++)
440  {
441  label nodeI;
442  faceLine >> nodeI;
443  f[2-i] = nodeToPoint[nodeI];
444  }
445 
446 
447  if (findFace(mesh, f) >= mesh.nInternalFaces())
448  {
449  boundaryFaces[faceI] = f;
450 
451  if (nFaceAttr > 0)
452  {
453  // First attribute is the region number
454  faceLine >> region;
455 
456 
457  // Get Foam patchID and update region->patch table.
458  label patchI = 0;
459 
460  Map<label>::iterator patchFind =
461  regionToPatch.find(region);
462 
463  if (patchFind == regionToPatch.end())
464  {
465  patchI = nPatches;
466 
467  Info<< "Mapping tetgen region " << region
468  << " to Foam patch "
469  << patchI << endl;
470 
471  regionToPatch.insert(region, nPatches++);
472  }
473  else
474  {
475  patchI = patchFind();
476  }
477 
478  boundaryPatch[faceI] = patchI;
479 
480  // Skip remaining attributes
481  for (label i = 1; i < nFaceAttr; i++)
482  {
483  faceLine >> dummy;
484  }
485  }
486 
487  faceI++;
488  }
489  }
490  }
491 
492 
493  // Trim
494  boundaryFaces.setSize(faceI);
495  boundaryPatch.setSize(faceI);
496 
497 
498  // Print region to patch mapping
499  Info<< "Regions:" << endl;
500 
501  for
502  (
503  Map<label>::const_iterator iter = regionToPatch.begin();
504  iter != regionToPatch.end();
505  ++iter
506  )
507  {
508  Info<< " region:" << iter.key() << '\t' << "patch:"
509  << iter() << endl;
510  }
511  Info<< endl;
512 
513 
514  // Storage for boundary faces
515  faceListList patchFaces(nPatches);
516  wordList patchNames(nPatches);
517 
518  forAll(patchNames, patchI)
519  {
520  patchNames[patchI] = word("patch") + name(patchI);
521  }
522 
523  wordList patchTypes(nPatches, polyPatch::typeName);
524  word defaultFacesName = "defaultFaces";
525  word defaultFacesType = polyPatch::typeName;
526  wordList patchPhysicalTypes(nPatches, polyPatch::typeName);
527 
528 
529  // Sort boundaryFaces by patch using boundaryPatch.
530  List<DynamicList<face> > allPatchFaces(nPatches);
531 
532  forAll(boundaryPatch, faceI)
533  {
534  label patchI = boundaryPatch[faceI];
535 
536  allPatchFaces[patchI].append(boundaryFaces[faceI]);
537  }
538 
539  Info<< "Patch sizes:" << endl;
540 
541  forAll(allPatchFaces, patchI)
542  {
543  Info<< " " << patchNames[patchI] << " : "
544  << allPatchFaces[patchI].size() << endl;
545 
546  patchFaces[patchI].transfer(allPatchFaces[patchI]);
547  }
548 
549  Info<< endl;
550 
551 
552  meshPtr.reset
553  (
554  new polyMesh
555  (
556  IOobject
557  (
559  runTime.constant(),
560  runTime
561  ),
562  xferMove(points),
563  cells,
564  patchFaces,
565  patchNames,
566  patchTypes,
570  )
571  );
572  }
573 
574 
575  Info<< "Writing mesh to " << runTime.constant() << endl << endl;
576 
577  meshPtr().write();
578 
579  Info<< "End\n" << endl;
580 
581  return 0;
582 }
583 
584 
585 // ************************ vim: set sw=4 sts=4 et: ************************ //