applications/utilities/mesh/conversion/netgenNeutralToFoam/netgenNeutralToFoam.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 neutral file format as written by Netgen v4.4.
  26
  27     Example:
  28
  29     9
  30       1.000000  1.000000  1.000000
  31       0.000000  1.000000  1.000000
  32       0.000000  0.000000  1.000000
  33       1.000000  0.000000  1.000000
  34       0.000000  1.000000  0.000000
  35       1.000000  1.000000  0.000000
  36       1.000000  0.000000  0.000000
  37       0.000000  0.000000  0.000000
  38       0.500000  0.500000  0.500000
  39     12
  40        1          7        8        9        3
  41        1          5        9        6        8
  42        1          5        9        2        1
  43        1          4        9        7        6
  44        1          7        8        6        9
  45        1          4        6        1        9
  46        1          5        9        8        2
  47        1          4        1        2        9
  48        1          1        6        5        9
  49        1          2        3        4        9
  50        1          8        9        3        2
  51        1          4        9        3        7
  52     12
  53        1            1        2        4
  54        1            3        4        2
  55        2            5        6        8
  56        2            7        8        6
  57        3            1        4        6
  58        3            7        6        4
  59        5            2        1        5
  60        5            6        5        1
  61        5            3        2        8
  62        5            5        8        2
  63        6            4        3        7
  64        6            8        7        3
  65
  66 NOTE:
  67     - reverse order of boundary faces using geometric test.
  68       (not very space efficient)
  69     - order of tet vertices only tested on one file.
  70     - all patch/cell/vertex numbers offset by one.
  71
  72 \*---------------------------------------------------------------------------*/
  73
  74 #include "argList.H"
  75 #include "Time.H"
  76 #include "polyMesh.H"
  77 #include "IFstream.H"
  78 #include "polyPatch.H"
  79 #include "cellModeller.H"
  80 #include "triFace.H"
  81
  82 using namespace Foam;
  83
  84 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  85
  86 // Main program:
  87
  88 int main(int argc, char *argv[])
  89 {
  90     argList::validArgs.append("Neutral file");
  91
  92 #   include "setRootCase.H"
  93 #   include "createTime.H"
  94
  95     IFstream str(args[1]);
  96
  97     //
  98     // Read nodes.
  99     //
 100     label nNodes(readLabel(str));
 101
 102     Info<< "nNodes:" << nNodes << endl;
 103
 104     pointField points(nNodes);
 105
 106     forAll(points, pointI)
 107     {
 108         scalar x,y,z;
 109
 110         str >> x >> y >> z;
 111
 112         points[pointI] = point(x, y, z);
 113     }
 114
 115
 116
 117
 118     label nTets(readLabel(str));
 119
 120     Info<< "nTets:" << nTets << endl;
 121
 122     const cellModel& tet = *(cellModeller::lookup("tet"));
 123
 124     cellShapeList cells(nTets);
 125
 126     labelList tetPoints(4);
 127
 128     forAll(cells, cellI)
 129     {
 130         label domain(readLabel(str));
 131
 132         if (domain != 1)
 133         {
 134             WarningIn(args.executable())
 135                 << "Cannot handle multiple domains"
 136                 << nl << "Ignoring domain " << domain << " setting on line "
 137                 << str.lineNumber() << endl;
 138         }
 139
 140         tetPoints[1] = readLabel(str) - 1;
 141         tetPoints[0] = readLabel(str) - 1;
 142         tetPoints[2] = readLabel(str) - 1;
 143         tetPoints[3] = readLabel(str) - 1;
 144
 145         cells[cellI] = cellShape(tet, tetPoints);
 146     }
 147
 148
 149
 150     label nFaces(readLabel(str));
 151
 152     Info<< "nFaces:" << nFaces << endl;
 153
 154     // Unsorted boundary faces
 155     faceList boundaryFaces(nFaces);
 156
 157     // For each boundary faces the Foam patchID
 158     labelList boundaryPatch(nFaces, -1);
 159
 160     // Max patch.
 161     label maxPatch = 0;
 162
 163     // Boundary faces as three vertices
 164     HashTable<label, triFace, Hash<triFace> > vertsToBoundary(nFaces);
 165
 166     forAll(boundaryFaces, faceI)
 167     {
 168         label patchI(readLabel(str));
 169
 170         if (patchI < 0)
 171         {
 172             FatalErrorIn(args.executable())
 173                 << "Invalid boundary region number " << patchI
 174                 << " on line " << str.lineNumber()
 175                 << exit(FatalError);
 176         }
 177
 178
 179         maxPatch = max(maxPatch, patchI);
 180
 181         triFace tri(readLabel(str)-1, readLabel(str)-1, readLabel(str)-1);
 182
 183         // Store boundary face as is for now. Later on reverse it.
 184         boundaryFaces[faceI].setSize(3);
 185         boundaryFaces[faceI][0] = tri[0];
 186         boundaryFaces[faceI][1] = tri[1];
 187         boundaryFaces[faceI][2] = tri[2];
 188         boundaryPatch[faceI] = patchI;
 189
 190         vertsToBoundary.insert(tri, faceI);
 191     }
 192
 193     label nPatches = maxPatch + 1;
 194
 195
 196     // Use hash of points to get owner cell and orient the boundary face.
 197     // For storage reasons I store the triangles and loop over the cells instead
 198     // of the other way around (store cells and loop over triangles) though
 199     // that would be faster.
 200     forAll(cells, cellI)
 201     {
 202         const cellShape& cll = cells[cellI];
 203
 204         // Get the four (outwards pointing) faces of the cell
 205         faceList tris(cll.faces());
 206
 207         forAll(tris, i)
 208         {
 209             const face& f = tris[i];
 210
 211             // Is there any boundary face with same vertices?
 212             // (uses commutative hash)
 213             HashTable<label, triFace, Hash<triFace> >::iterator iter =
 214                 vertsToBoundary.find(triFace(f[0], f[1], f[2]));
 215
 216             if (iter != vertsToBoundary.end())
 217             {
 218                 label faceI = iter();
 219                 const triFace& tri = iter.key();
 220
 221                 // Determine orientation of tri v.s. cell centre.
 222                 point cc(cll.centre(points));
 223                 point fc(tri.centre(points));
 224                 vector fn(tri.normal(points));
 225
 226                 if (((fc - cc) & fn) < 0)
 227                 {
 228                     // Boundary face points inwards. Flip.
 229                     boundaryFaces[faceI].flip();
 230                 }
 231
 232                 // Done this face so erase from hash
 233                 vertsToBoundary.erase(iter);
 234             }
 235         }
 236     }
 237
 238
 239     if (vertsToBoundary.size())
 240     {
 241         // Didn't find cells connected to boundary faces.
 242         WarningIn(args.executable())
 243             << "There are boundary faces without attached cells."
 244             << "Boundary faces (as triFaces):" << vertsToBoundary.toc()
 245             << endl;
 246     }
 247
 248
 249     // Storage for boundary faces sorted into patches
 250
 251     faceListList patchFaces(nPatches);
 252
 253     wordList patchNames(nPatches);
 254
 255     forAll(patchNames, patchI)
 256     {
 257         patchNames[patchI] = word("patch") + name(patchI);
 258     }
 259
 260     wordList patchTypes(nPatches, polyPatch::typeName);
 261     word defaultFacesName = "defaultFaces";
 262     word defaultFacesType = polyPatch::typeName;
 263     wordList patchPhysicalTypes(nPatches, polyPatch::typeName);
 264
 265     {
 266         // Sort boundaryFaces by patch.
 267         List<DynamicList<face> > allPatchFaces(nPatches);
 268
 269         forAll(boundaryPatch, faceI)
 270         {
 271             label patchI = boundaryPatch[faceI];
 272
 273             allPatchFaces[patchI].append(boundaryFaces[faceI]);
 274         }
 275
 276         Info<< "Patches:" << nl
 277             << "\tNeutral Boundary\tPatch name\tSize" << nl
 278             << "\t----------------\t----------\t----" << endl;
 279
 280         forAll(allPatchFaces, patchI)
 281         {
 282             Info<< '\t' << patchI << "\t\t\t"
 283                 << patchNames[patchI] << "\t\t"
 284                 << allPatchFaces[patchI].size() << endl;
 285
 286             patchFaces[patchI].transfer(allPatchFaces[patchI]);
 287         }
 288
 289         Info<< endl;
 290     }
 291
 292
 293     polyMesh mesh
 294     (
 295         IOobject
 296         (
 297             polyMesh::defaultRegion,
 298             runTime.constant(),
 299             runTime
 300         ),
 301         xferMove(points),
 302         cells,
 303         patchFaces,
 304         patchNames,
 305         patchTypes,
 306         defaultFacesName,
 307         defaultFacesType,
 308         patchPhysicalTypes
 309     );
 310
 311     Info<< "Writing mesh to " << runTime.constant() << endl << endl;
 312
 313     mesh.write();
 314
 315
 316     Info<< "End\n" << endl;
 317
 318     return 0;
 319 }
 320
 321
 322 // ************************************************************************* //