applications/utilities/mesh/conversion/tetgenToFoam/tetgenToFoam.C

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