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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2015-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 tetgenToFoam
29
30Group
31 grpMeshConversionUtilities
32
33Description
34 Convert tetgen .ele and .node and .face files to an OpenFOAM mesh.
35
36 Make sure to use add boundary attributes to the smesh file
37 (5 fifth column in the element section)
38 and run tetgen with -f option.
39
40 Sample smesh file:
41 \verbatim
42 # cube.smesh -- A 10x10x10 cube
43 8 3
44 1 0 0 0
45 2 0 10 0
46 3 10 10 0
47 4 10 0 0
48 5 0 0 10
49 6 0 10 10
50 7 10 10 10
51 8 10 0 10
52 6 1 # 1 for boundary info present
53 4 1 2 3 4 11 # region number 11
54 4 5 6 7 8 21 # region number 21
55 4 1 2 6 5 3
56 4 4 3 7 8 43
57 4 1 5 8 4 5
58 4 2 6 7 3 65
59 0
60 0
61 \endverbatim
62
63Note
64 - for some reason boundary faces point inwards. I just reverse them
65 always. Might use some geometric check instead.
66 - marked faces might not actually be boundary faces of mesh.
67 This is hopefully handled now by first constructing without boundaries
68 and then reconstructing with boundary faces.
69
70\*---------------------------------------------------------------------------*/
71
72#include "argList.H"
73#include "Time.H"
74#include "polyMesh.H"
75#include "IFstream.H"
76#include "cellModel.H"
77
78using namespace Foam;
79
80// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
81
82// Find label of face.
83label findFace(const primitiveMesh& mesh, const face& f)
84{
85 const labelList& pFaces = mesh.pointFaces()[f[0]];
86
87 forAll(pFaces, i)
88 {
89 label facei = pFaces[i];
90
91 if (mesh.faces()[facei] == f)
92 {
93 return facei;
94 }
95 }
96
98 << "Cannot find face " << f << " in mesh." << abort(FatalError);
99
100 return -1;
101}
102
103
104
105int main(int argc, char *argv[])
106{
107 argList::addNote
108 (
109 "Convert tetgen .ele and .node and .face files to an OpenFOAM mesh"
110 );
111
112 argList::addArgument("prefix", "The prefix for the input tetgen files");
113 argList::addBoolOption
114 (
115 "noFaceFile",
116 "Skip reading .face file for boundary information"
117 );
118
119 #include "setRootCase.H"
120 #include "createTime.H"
121
122 const auto prefix = args.get<fileName>(1);
123 const bool readFaceFile = !args.found("noFaceFile");
124
125 const fileName nodeFile(prefix + ".node");
126 const fileName eleFile(prefix + ".ele");
127 const fileName faceFile(prefix + ".face");
128
129 if (!readFaceFile)
130 {
131 Info<< "Files:" << endl
132 << " nodes : " << nodeFile << endl
133 << " elems : " << eleFile << endl
134 << endl;
135 }
136 else
137 {
138 Info<< "Files:" << endl
139 << " nodes : " << nodeFile << endl
140 << " elems : " << eleFile << endl
141 << " faces : " << faceFile << endl
142 << endl;
143
144 Info<< "Reading .face file for boundary information" << nl << endl;
145 }
146
147 if (!isFile(nodeFile) || !isFile(eleFile))
148 {
150 << "Cannot read " << nodeFile << " or " << eleFile
151 << exit(FatalError);
152 }
153
154 if (readFaceFile && !isFile(faceFile))
155 {
157 << "Cannot read " << faceFile << endl
158 << "Did you run tetgen with -f option?" << endl
159 << "If you don't want to read the .face file and thus not have"
160 << " patches please\nrerun with the -noFaceFile option"
161 << exit(FatalError);
162 }
163
164
165 IFstream nodeStream(nodeFile);
166
167 //
168 // Read nodes.
169 //
170
171 // Read header.
172 string line;
173
174 do
175 {
176 nodeStream.getLine(line);
177 }
178 while (line.starts_with('#'));
179
180 IStringStream nodeLine(line);
181
182 label nNodes, nDims, nNodeAttr;
183 bool hasRegion;
184
185 nodeLine >> nNodes >> nDims >> nNodeAttr >> hasRegion;
186
187
188 Info<< "Read .node header:" << endl
189 << " nodes : " << nNodes << endl
190 << " nDims : " << nDims << endl
191 << " nAttr : " << nNodeAttr << endl
192 << " hasRegion : " << hasRegion << endl
193 << endl;
194
195 //
196 // read points
197 //
198
199 pointField points(nNodes);
200 Map<label> nodeToPoint(nNodes);
201
202 {
203 labelList pointIndex(nNodes);
204
205 label pointi = 0;
206
207 while (nodeStream.good())
208 {
209 nodeStream.getLine(line);
210
211 if (line.size() && line[0] != '#')
212 {
213 IStringStream nodeLine(line);
214
215 label nodeI;
216 scalar x, y, z;
217 label dummy;
218
219 nodeLine >> nodeI >> x >> y >> z;
220
221 for (label i = 0; i < nNodeAttr; i++)
222 {
223 nodeLine >> dummy;
224 }
225
226 if (hasRegion)
227 {
228 nodeLine >> dummy;
229 }
230
231 // Store point and node number.
232 points[pointi] = point(x, y, z);
233 nodeToPoint.insert(nodeI, pointi);
234 pointi++;
235 }
236 }
237 if (pointi != nNodes)
238 {
239 FatalIOErrorInFunction(nodeStream)
240 << "Only " << pointi << " nodes present instead of " << nNodes
241 << " from header." << exit(FatalIOError);
242 }
243 }
244
245 //
246 // read elements
247 //
248
249 IFstream eleStream(eleFile);
250
251 do
252 {
253 eleStream.getLine(line);
254 }
255 while (line.starts_with('#'));
256
257 IStringStream eleLine(line);
258
259 label nTets, nPtsPerTet, nElemAttr;
260
261 eleLine >> nTets >> nPtsPerTet >> nElemAttr;
262
263
264 Info<< "Read .ele header:" << endl
265 << " tets : " << nTets << endl
266 << " pointsPerTet : " << nPtsPerTet << endl
267 << " nAttr : " << nElemAttr << endl
268 << endl;
269
270 if (nPtsPerTet != 4)
271 {
272 FatalIOErrorInFunction(eleStream)
273 << "Cannot handle tets with "
274 << nPtsPerTet << " points per tetrahedron in .ele file" << endl
275 << "Can only handle tetrahedra with four points"
276 << exit(FatalIOError);
277 }
278
279 if (nElemAttr != 0)
280 {
282 << "Element attributes (third elemenent in .ele header)"
283 << " not used" << endl;
284 }
285
286
287 const cellModel& tet = cellModel::ref(cellModel::TET);
288
290
291 cellShapeList cells(nTets);
292 label celli = 0;
293
294 while (eleStream.good())
295 {
296 eleStream.getLine(line);
297
298 if (line.size() && line[0] != '#')
299 {
300 IStringStream eleLine(line);
301
302 label elemI;
303 eleLine >> elemI;
304
305 for (label i = 0; i < 4; i++)
306 {
307 label nodeI;
308 eleLine >> nodeI;
309 tetPoints[i] = nodeToPoint[nodeI];
310 }
311
312 cells[celli++].reset(tet, tetPoints);
313
314 // Skip attributes
315 for (label i = 0; i < nElemAttr; i++)
316 {
317 label dummy;
318
319 eleLine >> dummy;
320 }
321 }
322 }
323
324
325 //
326 // Construct mesh with default boundary only
327 //
328
330 (
332 (
333 polyMesh::defaultRegion,
334 runTime.constant(),
335 runTime,
336 IOobject::NO_READ,
337 IOobject::AUTO_WRITE
338 ),
339 pointField(points), // Copy of points
340 cells,
341 faceListList(),
342 wordList(), // boundaryPatchNames
343 wordList(), // boundaryPatchTypes
344 "defaultFaces",
345 polyPatch::typeName,
346 wordList()
347 );
348 const polyMesh& mesh = *meshPtr;
349
350
351 if (readFaceFile)
352 {
353 label nPatches = 0;
354
355 // List of Foam vertices per boundary face
356 faceList boundaryFaces;
357
358 // For each boundary faces the Foam patchID
360
361 //
362 // read boundary faces
363 //
364
365 IFstream faceStream(faceFile);
366
367 do
368 {
369 faceStream.getLine(line);
370 }
371 while (line.starts_with('#'));
372
373 IStringStream faceLine(line);
374
375 label nFaces, nFaceAttr;
376
377 faceLine >> nFaces >> nFaceAttr;
378
379
380 Info<< "Read .face header:" << endl
381 << " faces : " << nFaces << endl
382 << " nAttr : " << nFaceAttr << endl
383 << endl;
384
385
386 if (nFaceAttr != 1)
387 {
388 FatalIOErrorInFunction(faceStream)
389 << "Expect boundary markers to be"
390 << " present in .face file." << endl
391 << "This is the second number in the header which is now:"
392 << nFaceAttr << exit(FatalIOError);
393 }
394
395 // List of Foam vertices per boundary face
396 boundaryFaces.setSize(nFaces);
397
398 // For each boundary faces the Foam patchID
399 boundaryPatch.setSize(nFaces);
400 boundaryPatch = -1;
401
402 label facei = 0;
403
404 // Region to patch conversion
405 Map<label> regionToPatch;
406
407 face f(3);
408
409 while (faceStream.good())
410 {
411 faceStream.getLine(line);
412
413 if (line.size() && line[0] != '#')
414 {
415 IStringStream faceLine(line);
416
417 label tetGenFacei, dummy, region;
418
419 faceLine >> tetGenFacei;
420
421 // Read face and reverse orientation (Foam needs outwards
422 // pointing)
423 for (label i = 0; i < 3; i++)
424 {
425 label nodeI;
426 faceLine >> nodeI;
427 f[2-i] = nodeToPoint[nodeI];
428 }
429
430
431 if (findFace(mesh, f) >= mesh.nInternalFaces())
432 {
433 boundaryFaces[facei] = f;
434
435 if (nFaceAttr > 0)
436 {
437 // First attribute is the region number
438 faceLine >> region;
439
440
441 // Get Foam patchID and update region->patch table.
442 label patchi = 0;
443
444 const auto patchFind = regionToPatch.cfind(region);
445
446 if (patchFind.found())
447 {
448 patchi = *patchFind;
449 }
450 else
451 {
452 patchi = nPatches;
453
454 Info<< "Mapping tetgen region " << region
455 << " to patch "
456 << patchi << endl;
457
458 regionToPatch.insert(region, nPatches++);
459 }
460
461 boundaryPatch[facei] = patchi;
462
463 // Skip remaining attributes
464 for (label i = 1; i < nFaceAttr; i++)
465 {
466 faceLine >> dummy;
467 }
468 }
469
470 facei++;
471 }
472 }
473 }
474
475
476 // Trim
477 boundaryFaces.setSize(facei);
478 boundaryPatch.setSize(facei);
479
480
481 // Print region to patch mapping
482 Info<< "Regions:" << endl;
483
484 forAllConstIters(regionToPatch, iter)
485 {
486 Info<< " region:" << iter.key()
487 << '\t' << "patch:" << iter.val() << nl;
488 }
489 Info<< endl;
490
491
492 // Storage for boundary faces
493 faceListList patchFaces(nPatches);
495
496 forAll(patchNames, patchi)
497 {
498 patchNames[patchi] = polyPatch::defaultName(patchi);
499 }
500
501 wordList patchTypes(nPatches, polyPatch::typeName);
502 word defaultFacesName = "defaultFaces";
503 word defaultFacesType = polyPatch::typeName;
504 wordList patchPhysicalTypes(nPatches, polyPatch::typeName);
505
506
507 // Sort boundaryFaces by patch using boundaryPatch.
508 List<DynamicList<face>> allPatchFaces(nPatches);
509
510 forAll(boundaryPatch, facei)
511 {
512 label patchi = boundaryPatch[facei];
513
514 allPatchFaces[patchi].append(boundaryFaces[facei]);
515 }
516
517 Info<< "Patch sizes:" << endl;
518
519 forAll(allPatchFaces, patchi)
520 {
521 Info<< " " << patchNames[patchi] << " : "
522 << allPatchFaces[patchi].size() << endl;
523
524 patchFaces[patchi].transfer(allPatchFaces[patchi]);
525 }
526
527 Info<< endl;
528
529
530 meshPtr.reset
531 (
532 new polyMesh
533 (
535 (
536 polyMesh::defaultRegion,
537 runTime.constant(),
538 runTime
539 ),
540 std::move(points),
541 cells,
542 patchFaces,
547 patchPhysicalTypes
548 )
549 );
550 }
551
552 // Set the precision of the points data to 10
553 IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
554
555 Info<< "Writing mesh to " << runTime.constant() << endl << endl;
556
557 meshPtr().removeFiles();
558 meshPtr().write();
559
560 Info<< "End\n" << endl;
561
562 return 0;
563}
564
565
566// ************************************************************************* //
scalar y
Y[inertIndex] max(0.0)
const_iterator cfind(const Key &key) const
Find and return an const_iterator set at the hashed entry.
Definition: HashTableI.H:141
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:180
Input from file stream, using an ISstream.
Definition: IFstream.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
Input from string buffer, using a ISstream. Always UNCOMPRESSED.
Definition: StringStream.H:112
void setSize(const label n)
Alias for resize()
Definition: List.H:218
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
Like polyPatch but without reference to mesh. Used in boundaryMesh to hold data on patches....
Definition: boundaryPatch.H:57
Maps a geometry to a set of cell primitives.
Definition: cellModel.H:73
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A class for handling file names.
Definition: fileName.H:76
A line primitive.
Definition: line.H:68
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:79
Tet storage. Null constructable (unfortunately tetrahedron<point, point> is not)
Definition: tetPoints.H:56
A class for handling words, derived from Foam::string.
Definition: word.H:68
dynamicFvMesh & mesh
engineTime & runTime
Foam::autoPtr< Foam::dynamicFvMesh > meshPtr
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:473
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
const label nPatches
const pointField & points
const cellShapeList & cells
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
errorManip< error > abort(error &err)
Definition: errorManip.H:144
bool isFile(const fileName &name, const bool checkGzip=true, const bool followLink=true)
Does the name exist as a FILE in the file system?
Definition: MSwindows.C:666
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
wordList patchTypes(nPatches)
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
wordList patchNames(nPatches)
labelList f(nPoints)
word defaultFacesName
Definition: readKivaGrid.H:455
word defaultFacesType
Definition: readKivaGrid.H:456
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278