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