netgenNeutralToFoam.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.
 
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.
 
    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
Application
    netgenNeutralToFoam
 
Description
    Converts neutral file format as written by Netgen v4.4.
 
    Example:
 
    9
      1.000000  1.000000  1.000000
      0.000000  1.000000  1.000000
      0.000000  0.000000  1.000000
      1.000000  0.000000  1.000000
      0.000000  1.000000  0.000000
      1.000000  1.000000  0.000000
      1.000000  0.000000  0.000000
      0.000000  0.000000  0.000000
      0.500000  0.500000  0.500000
    12
       1          7        8        9        3
       1          5        9        6        8
       1          5        9        2        1
       1          4        9        7        6
       1          7        8        6        9
       1          4        6        1        9
       1          5        9        8        2
       1          4        1        2        9
       1          1        6        5        9
       1          2        3        4        9
       1          8        9        3        2
       1          4        9        3        7
    12
       1            1        2        4
       1            3        4        2
       2            5        6        8
       2            7        8        6
       3            1        4        6
       3            7        6        4
       5            2        1        5
       5            6        5        1
       5            3        2        8
       5            5        8        2
       6            4        3        7
       6            8        7        3
 
NOTE:
    - reverse order of boundary faces using geometric test.
      (not very space efficient)
    - order of tet vertices only tested on one file.
    - all patch/cell/vertex numbers offset by one.
 
\*---------------------------------------------------------------------------*/
 
#include "argList.H"
#include "Time.H"
#include "polyMesh.H"
#include "IFstream.H"
#include "polyPatch.H"
#include "cellModeller.H"
#include "triFace.H"
 
using namespace Foam;
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 
int main(int argc, char *argv[])
{
    argList::validArgs.append("Neutral file");
 
    #include "setRootCase.H"
    #include "createTime.H"
 
    IFstream str(args[1]);
 
    //
    // Read nodes.
    //
    label nNodes(readLabel(str));
 
    Info<< "nNodes:" << nNodes << endl;
 
    pointField points(nNodes);
 
    forAll(points, pointi)
    {
        scalar x,y,z;
 
        str >> x >> y >> z;
 
        points[pointi] = point(x, y, z);
    }
 
 
 
 
    label nTets(readLabel(str));
 
    Info<< "nTets:" << nTets << endl;
 
    const cellModel& tet = *(cellModeller::lookup("tet"));
 
    cellShapeList cells(nTets);
 
    labelList tetPoints(4);
 
    forAll(cells, celli)
    {
        label domain(readLabel(str));
 
        if (domain != 1)
        {
            WarningInFunction
                << "Cannot handle multiple domains"
                << nl << "Ignoring domain " << domain << " setting on line "
                << str.lineNumber() << endl;
        }
 
        tetPoints[1] = readLabel(str) - 1;
        tetPoints[0] = readLabel(str) - 1;
        tetPoints[2] = readLabel(str) - 1;
        tetPoints[3] = readLabel(str) - 1;
 
        cells[celli] = cellShape(tet, tetPoints);
    }
 
 
 
    label nFaces(readLabel(str));
 
    Info<< "nFaces:" << nFaces << endl;
 
    // Unsorted boundary faces
    faceList boundaryFaces(nFaces);
 
    // For each boundary faces the Foam patchID
    labelList boundaryPatch(nFaces, -1);
 
    // Max patch.
    label maxPatch = 0;
 
    // Boundary faces as three vertices
    HashTable<label, triFace, Hash<triFace>> vertsToBoundary(nFaces);
 
    forAll(boundaryFaces, facei)
    {
        label patchi(readLabel(str));
 
        if (patchi < 0)
        {
            FatalErrorInFunction
                << "Invalid boundary region number " << patchi
                << " on line " << str.lineNumber()
                << exit(FatalError);
        }
 
 
        maxPatch = max(maxPatch, patchi);
 
        triFace tri(readLabel(str)-1, readLabel(str)-1, readLabel(str)-1);
 
        // Store boundary face as is for now. Later on reverse it.
        boundaryFaces[facei].setSize(3);
        boundaryFaces[facei][0] = tri[0];
        boundaryFaces[facei][1] = tri[1];
        boundaryFaces[facei][2] = tri[2];
        boundaryPatch[facei] = patchi;
 
        vertsToBoundary.insert(tri, facei);
    }
 
    label nPatches = maxPatch + 1;
 
 
    // Use hash of points to get owner cell and orient the boundary face.
    // For storage reasons I store the triangles and loop over the cells instead
    // of the other way around (store cells and loop over triangles) though
    // that would be faster.
    forAll(cells, celli)
    {
        const cellShape& cll = cells[celli];
 
        // Get the four (outwards pointing) faces of the cell
        faceList tris(cll.faces());
 
        forAll(tris, i)
        {
            const face& f = tris[i];
 
            // Is there any boundary face with same vertices?
            // (uses commutative hash)
            HashTable<label, triFace, Hash<triFace>>::iterator iter =
                vertsToBoundary.find(triFace(f[0], f[1], f[2]));
 
            if (iter != vertsToBoundary.end())
            {
                label facei = iter();
                const triFace& tri = iter.key();
 
                // Determine orientation of tri v.s. cell centre.
                point cc(cll.centre(points));
                point fc(tri.centre(points));
                vector fn(tri.area(points));
 
                if (((fc - cc) & fn) < 0)
                {
                    // Boundary face points inwards. Flip.
                    boundaryFaces[facei].flip();
                }
 
                // Done this face so erase from hash
                vertsToBoundary.erase(iter);
            }
        }
    }
 
 
    if (vertsToBoundary.size())
    {
        // Didn't find cells connected to boundary faces.
        WarningInFunction
            << "There are boundary faces without attached cells."
            << "Boundary faces (as triFaces):" << vertsToBoundary.toc()
            << endl;
    }
 
 
    // Storage for boundary faces sorted into patches
 
    faceListList patchFaces(nPatches);
 
    wordList patchNames(nPatches);
 
    forAll(patchNames, patchi)
    {
        patchNames[patchi] = word("patch") + name(patchi);
    }
 
    wordList patchTypes(nPatches, polyPatch::typeName);
    word defaultFacesName = "defaultFaces";
    word defaultFacesType = polyPatch::typeName;
    wordList patchPhysicalTypes(nPatches, polyPatch::typeName);
 
    {
        // Sort boundaryFaces by patch.
        List<DynamicList<face>> allPatchFaces(nPatches);
 
        forAll(boundaryPatch, facei)
        {
            label patchi = boundaryPatch[facei];
 
            allPatchFaces[patchi].append(boundaryFaces[facei]);
        }
 
        Info<< "Patches:" << nl
            << "\tNeutral Boundary\tPatch name\tSize" << nl
            << "\t----------------\t----------\t----" << endl;
 
        forAll(allPatchFaces, patchi)
        {
            Info<< '\t' << patchi << "\t\t\t"
                << patchNames[patchi] << "\t\t"
                << allPatchFaces[patchi].size() << endl;
 
            patchFaces[patchi].transfer(allPatchFaces[patchi]);
        }
 
        Info<< endl;
    }
 
 
    polyMesh mesh
    (
        IOobject
        (
            polyMesh::defaultRegion,
            runTime.constant(),
            runTime
        ),
        move(points),
        cells,
        patchFaces,
        patchNames,
        patchTypes,
        defaultFacesName,
        defaultFacesType,
        patchPhysicalTypes
    );
 
    Info<< "Writing mesh to " << runTime.constant() << endl << endl;
 
    mesh.write();
 
 
    Info<< "End\n" << endl;
 
    return 0;
}
 
 
// ************************************************************************* //