STARCDMeshReader.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2016-2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
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/>.

\*---------------------------------------------------------------------------*/

#include "STARCDMeshReader.H"
#include "oldCyclicPolyPatch.H"
#include "emptyPolyPatch.H"
#include "wallPolyPatch.H"
#include "symmetryPolyPatch.H"
#include "cellModel.H"
#include "ListOps.H"
#include "stringOps.H"
#include "IFstream.H"
#include "IOMap.H"

// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //

namespace Foam
{

    // Read and discard to newline
    static inline void readToNewline(ISstream& is)
    {
        char ch = '\n';
        do
        {
            is.get(ch);
        }
        while ((is) && ch != '\n');
    }

} // End namespace Foam


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

void Foam::fileFormats::STARCDMeshReader::readAux
(
    const objectRegistry& registry
)
{
    boundaryRegion_.readDict(registry);
    cellTable_.readDict(registry);
}


// Read points from <.vrt> file
//
/*---------------------------------------------------------------------------*\
Line 1:
  PROSTAR_VERTEX [newline]

Line 2:
  <version> 0 0 0 0 0 0 0 [newline]

Body:
  <vertexId>  <x>  <y>  <z> [newline]

\*---------------------------------------------------------------------------*/
Foam::label Foam::fileFormats::STARCDMeshReader::readPoints
(
    const fileName& inputName,
    const scalar scaleFactor
)
{
    label nPoints = 0, maxId = 0;
    token tok;

    // Pass 1:
    // get # points and maximum vertex label
    {
        IFstream is(inputName);
        readHeader(is, STARCDCore::HEADER_VRT);

        scalar x, y, z;

        while (is.read(tok).good() && tok.isLabel())
        {
            const label starVertexId = tok.labelToken();

            is >> x >> y >> z;

            maxId = max(maxId, starVertexId);
            ++nPoints;
        }
    }

    if (!nPoints)
    {
        FatalErrorInFunction
            << "No points in file " << inputName << nl
            << abort(FatalError);
    }

    Info<< "Number of points  = " << nPoints << endl;

    // Set sizes and reset to invalid values

    points_.setSize(nPoints);
    mapToFoamPointId_.setSize(maxId+1);

    //- Original Point number for a given vertex
    // might need again in the future

    mapToFoamPointId_ = -1;

    // Pass 2:
    // construct pointList and conversion table
    // from Star vertex numbers to Foam point labels
    {
        IFstream is(inputName);
        readHeader(is, STARCDCore::HEADER_VRT);

        label pointi = 0;
        while (is.read(tok).good() && tok.isLabel())
        {
            const label starVertexId = tok.labelToken();

            is  >> points_[pointi].x()
                >> points_[pointi].y()
                >> points_[pointi].z();

            // might need again in the future
            mapToFoamPointId_[starVertexId] = pointi;
            ++pointi;
        }

        if (nPoints > pointi)
        {
            nPoints = pointi;
            points_.setSize(nPoints);
            // might need again in the future
        }

        if
        (
            scaleFactor > 0
         && (scaleFactor > 1.0 + SMALL || scaleFactor < 1.0 - SMALL)
        )
        {
            points_ *= scaleFactor;
        }
    }

    return maxId;
}


// Read cells from <.cel> file
//
/*---------------------------------------------------------------------------*\
Line 1:
  PROSTAR_CELL [newline]

Line 2:
  <version> 0 0 0 0 0 0 0 [newline]

Body:
  <cellId>  <shapeId>  <nLabels>  <cellTableId>  <typeId> [newline]
  <cellId>  <int1> .. <int8>
  <cellId>  <int9> .. <int16>

 with shapeId:
 *   1 = point
 *   2 = line
 *   3 = shell
 *  11 = hexa
 *  12 = prism
 *  13 = tetra
 *  14 = pyramid
 * 255 = polyhedron

 with typeId
 *   1 = fluid
 *   2 = solid
 *   3 = baffle
 *   4 = shell
 *   5 = line
 *   6 = point

For primitive cell shapes, the number of vertices will never exceed 8 (hexa)
and corresponds to <nLabels>.
For polyhedral, <nLabels> includes an index table comprising beg/end pairs
for each cell face.

Strictly speaking, we only need the cellModeller for adding boundaries.
\*---------------------------------------------------------------------------*/

void Foam::fileFormats::STARCDMeshReader::readCells(const fileName& inputName)
{
    label nFluids = 0, nSolids = 0, nBaffles = 0, nShells = 0;
    label maxId = 0;
    token tok;

    bool unknownVertices = false;


    // Pass 1:
    // count nFluids, nSolids, nBaffle, nShell and maxId
    // also see if polyhedral cells were used
    {
        IFstream is(inputName);
        readHeader(is, STARCDCore::HEADER_CEL);

        label shapeId, nLabels, cellTableId, typeId;

        while (is.read(tok).good() && tok.isLabel())
        {
            const label starCellId = tok.labelToken();

            is  >> shapeId
                >> nLabels
                >> cellTableId
                >> typeId;

            // Skip the rest of the line
            readToNewline(is);

            // Max 8 indices per line
            while (nLabels > 0)
            {
                readToNewline(is);
                nLabels -= 8;
            }

            if (typeId == STARCDCore::starcdFluidType)
            {
                ++nFluids;
                maxId = max(maxId, starCellId);

                if (!cellTable_.found(cellTableId))
                {
                    cellTable_.setName(cellTableId);
                    cellTable_.setMaterial(cellTableId, "fluid");
                }
            }
            else if (typeId == STARCDCore::starcdSolidType)
            {
                ++nSolids;
                if (keepSolids_)
                {
                    maxId = max(maxId, starCellId);
                }

                if (!cellTable_.found(cellTableId))
                {
                    cellTable_.setName(cellTableId);
                    cellTable_.setMaterial(cellTableId, "solid");
                }
            }
            else if (typeId == STARCDCore::starcdBaffleType)
            {
                // baffles have no cellTable entry
                ++nBaffles;
                maxId = max(maxId, starCellId);
            }
            else if (typeId == STARCDCore::starcdShellType)
            {
                ++nShells;
                if (!cellTable_.found(cellTableId))
                {
                    cellTable_.setName(cellTableId);
                    cellTable_.setMaterial(cellTableId, "shell");
                }
            }
        }
    }

    const label nCells = nFluids + (keepSolids_ ? nSolids : 0);

    Info<< "Number of fluids  = " << nFluids << nl
        << "Number of baffles = " << nBaffles << nl
        << "Number of solids  = " << nSolids
        << (keepSolids_ ? " (treat as fluid)" : " (ignored)") << nl
        << "Number of shells  = " << nShells << " (ignored)" << nl;

    if (!nCells)
    {
        OSstream& err = FatalErrorInFunction;

        err << "No cells in file " << inputName << nl;

        if (nShells)
        {
            err << "Consists of shells only (typeId=4)." << nl;
        }

        err << nl
            << abort(FatalError);
    }

    cellFaces_.setSize(nCells);
    cellShapes_.setSize(nCells);
    cellTableId_.setSize(nCells);

    // information for the interfaces
    baffleFaces_.setSize(nBaffles);

    // extra space for baffles
    origCellId_.setSize(nCells + nBaffles);
    mapToFoamCellId_.setSize(maxId+1);
    mapToFoamCellId_ = -1;


    // avoid undefined shapes for polyhedra
    cellShape genericShape
    (
        cellModel::ref(cellModel::UNKNOWN), labelList()
    );

    // Pass 2:
    // construct cellFaces_ and possibly cellShapes_
    {
        IFstream is(inputName);
        readHeader(is, STARCDCore::HEADER_CEL);

        labelList starLabels(64);
        label ignoredLabel, shapeId, nLabels, cellTableId, typeId;

        label celli = 0, bafflei = 0;

        while (is.read(tok).good() && tok.isLabel())
        {
            const label starCellId = tok.labelToken();

            is  >> shapeId
                >> nLabels
                >> cellTableId
                >> typeId;

            if (nLabels > starLabels.size())
            {
                starLabels.setSize(nLabels);
            }
            starLabels = -1;

            // Read indices - max 8 per line
            for (label i = 0; i < nLabels; ++i)
            {
                if ((i % 8) == 0)
                {
                    is >> ignoredLabel; // Skip cellId for continuation lines
                }
                is >> starLabels[i];
            }

            // Skip solid cells
            if
            (
                typeId == STARCDCore::starcdSolidType
             && !keepSolids_
            )
            {
                continue;
            }

            // Determine the OpenFOAM cell shape
            const cellModel* curModelPtr = nullptr;

            // fluid/solid cells
            switch (shapeId)
            {
                case STARCDCore::starcdHex:
                    curModelPtr = cellModel::ptr(cellModel::HEX);
                    break;
                case STARCDCore::starcdPrism:
                    curModelPtr = cellModel::ptr(cellModel::PRISM);
                    break;
                case STARCDCore::starcdTet:
                    curModelPtr = cellModel::ptr(cellModel::TET);
                    break;
                case STARCDCore::starcdPyr:
                    curModelPtr = cellModel::ptr(cellModel::PYR);
                    break;
            }

            if (curModelPtr)
            {
                // primitive cell - use shapes

                // convert orig vertex Id to point label
                bool isBad = false;
                for (label i=0; i < nLabels; ++i)
                {
                    label pointId = mapToFoamPointId_[starLabels[i]];
                    if (pointId < 0)
                    {
                        Info<< "Cells inconsistent with vertex file. "
                            << "Star vertex " << starLabels[i]
                            << " does not exist" << endl;
                        isBad = true;
                        unknownVertices = true;
                    }
                    starLabels[i] = pointId;
                }

                if (isBad)
                {
                    continue;
                }

                // record original cell number and lookup
                origCellId_[celli] = starCellId;
                mapToFoamCellId_[starCellId] = celli;

                cellTableId_[celli] = cellTableId;
                cellShapes_[celli] = cellShape
                (
                    *curModelPtr,
                    SubList<label>(starLabels, nLabels)
                );

                cellFaces_[celli] = cellShapes_[celli].faces();
                ++celli;
            }
            else if (shapeId == STARCDCore::starcdPoly)
            {
                // polyhedral cell
                label nFaces = starLabels[0] - 1;

                // convert orig vertex id to point label
                // start with offset (skip the index table)
                bool isBad = false;
                for (label i=starLabels[0]; i < nLabels; ++i)
                {
                    label pointId = mapToFoamPointId_[starLabels[i]];
                    if (pointId < 0)
                    {
                        Info<< "Cells inconsistent with vertex file. "
                            << "Star vertex " << starLabels[i]
                            << " does not exist" << endl;
                        isBad = true;
                        unknownVertices = true;
                    }
                    starLabels[i] = pointId;
                }

                if (isBad)
                {
                    continue;
                }

                // traverse beg/end indices
                faceList faces(nFaces);
                label facei = 0;
                for (label i=0; i < nFaces; ++i)
                {
                    label beg = starLabels[i];
                    label n   = starLabels[i+1] - beg;

                    face f
                    (
                        SubList<label>(starLabels, n, beg)
                    );

                    f.collapse();

                    // valid faces only
                    if (f.size() >= 3)
                    {
                        faces[facei++] = f;
                    }
                }

                if (nFaces > facei)
                {
                    Info<< "star cell " << starCellId << " has "
                        << (nFaces - facei)
                        << " empty faces - could cause boundary "
                        << "addressing problems"
                        << endl;

                    nFaces = facei;
                    faces.setSize(nFaces);
                }

                if (nFaces < 4)
                {
                    FatalErrorInFunction
                        << "star cell " << starCellId << " has " << nFaces
                        << abort(FatalError);
                }

                // record original cell number and lookup
                origCellId_[celli] = starCellId;
                mapToFoamCellId_[starCellId] = celli;

                cellTableId_[celli] = cellTableId;
                cellShapes_[celli]  = genericShape;
                cellFaces_[celli]   = faces;
                ++celli;
            }
            else if (typeId == STARCDCore::starcdBaffleType)
            {
                // baffles

                // convert orig vertex id to point label
                bool isBad = false;
                for (label i=0; i < nLabels; ++i)
                {
                    label pointId = mapToFoamPointId_[starLabels[i]];
                    if (pointId < 0)
                    {
                        Info<< "Baffles inconsistent with vertex file. "
                            << "Star vertex " << starLabels[i]
                            << " does not exist" << endl;
                        isBad = true;
                        unknownVertices = true;
                    }
                    starLabels[i] = pointId;
                }

                if (isBad)
                {
                    continue;
                }


                face f
                (
                    SubList<label>(starLabels, nLabels)
                );

                f.collapse();

                // valid faces only
                if (f.size() >= 3)
                {
                    baffleFaces_[bafflei] = f;
                    // insert lookup addressing in normal list
                    mapToFoamCellId_[starCellId]  = nCells + bafflei;
                    origCellId_[nCells + bafflei] = starCellId;
                    ++bafflei;
                }
            }
        }

        baffleFaces_.setSize(bafflei);
    }

    if (unknownVertices)
    {
        FatalErrorInFunction
            << "cells with unknown vertices"
            << abort(FatalError);
    }

    // truncate lists

#ifdef DEBUG_READING
    Info<< "CELLS READ" << endl;
#endif

    // cleanup
    mapToFoamPointId_.clear();
}


// Read boundaries from <.bnd> file
//
/*---------------------------------------------------------------------------*\
Line 1:
  PROSTAR_BOUNDARY [newline]

Line 2:
  <version> 0 0 0 0 0 0 0 [newline]

Body:
  <boundId>  <cellId>  <cellFace>  <regionId>  0  <boundaryType> [newline]

where boundaryType is truncated to 4 characters from one of the following:
INLET
PRESSSURE
OUTLET
BAFFLE
etc,
\*---------------------------------------------------------------------------*/

void Foam::fileFormats::STARCDMeshReader::readBoundary
(
    const fileName& inputName
)
{
    label nPatches = 0, nFaces = 0, nBafflePatches = 0, maxId = 0;
    label starCellId, cellFaceId, starRegion, configNumber;
    token tok;
    word patchType;

    labelList mapToFoamPatchId(1000, label(-1));
    labelList nPatchFaces(1000, Zero);
    labelList origRegion(1000, Zero);
    patchTypes_.setSize(1000);

    //
    // Mapping between OpenFOAM and PROSTAR primitives
    // - needed for face mapping
    //
    const Map<label> shapeLookup =
    {
        { cellModel::ref(cellModel::HEX).index(), STARCDCore::starcdHex },
        { cellModel::ref(cellModel::PRISM).index(), STARCDCore::starcdPrism },
        { cellModel::ref(cellModel::TET).index(), STARCDCore::starcdTet },
        { cellModel::ref(cellModel::PYR).index(), STARCDCore::starcdPyr },
    };

    // Pass 1:
    // collect
    // no. of faces (nFaces), no. of patches (nPatches)
    // and for each of these patches the number of faces
    // (nPatchFaces[patchLabel])
    //
    // and a conversion table from Star regions to (Foam) patchLabels
    //
    // additionally note the no. of baffle patches (nBafflePatches)
    // so that we sort these to the end of the patch list
    // - this makes it easier to transfer them to an adjacent patch if reqd
    {
        IFstream is(inputName);

        if (is.good())
        {
            readHeader(is, STARCDCore::HEADER_BND);

            while (is.read(tok).good() && tok.isLabel())
            {
                // Ignore boundary id (not needed)

                ++nFaces;

                is  >> starCellId
                    >> cellFaceId
                    >> starRegion
                    >> configNumber
                    >> patchType;

                // Build translation table to convert star patch to foam patch
                label patchLabel = mapToFoamPatchId[starRegion];
                if (patchLabel == -1)
                {
                    patchLabel = nPatches;
                    mapToFoamPatchId[starRegion] = patchLabel;
                    origRegion[patchLabel] = starRegion;
                    patchTypes_[patchLabel] = patchType;

                    maxId = max(maxId, starRegion);

                    // should actually be case-insensitive
                    if (patchType == "BAFF")
                    {
                        ++nBafflePatches;
                    }
                    ++nPatches;
                }

                ++nPatchFaces[patchLabel];
            }

            if (nPatches == 0)
            {
                Info<< "No boundary faces in file " << inputName << endl;
            }
        }
        else
        {
            Info<< "Could not read boundary file " << inputName << endl;
        }
    }

    // keep empty patch region in reserve
    ++nPatches;
    Info<< "Number of patches = " << nPatches
        << " (including extra for missing) with "
        << nFaces << " faces" << endl;

    // resize
    origRegion.setSize(nPatches);
    patchTypes_.setSize(nPatches);
    patchNames_.setSize(nPatches);
    nPatchFaces.setSize(nPatches);

    // add our empty patch
    origRegion[nPatches-1] = 0;
    nPatchFaces[nPatches-1] = 0;
    patchTypes_[nPatches-1] = "none";

    // create names
    // - use 'Label' entry from "constant/boundaryRegion" dictionary
    forAll(patchTypes_, patchi)
    {
        bool fndName = false, fndType = false;

        auto iter = boundaryRegion_.cfind(origRegion[patchi]);

        if (iter.good())
        {
            const dictionary& dict = *iter;

            fndType = dict.readIfPresent("BoundaryType", patchTypes_[patchi]);
            fndName = dict.readIfPresent("Label", patchNames_[patchi]);
        }

        // Consistent names. Long form and in lowercase
        if (!fndType)
        {
            stringOps::inplaceLower(patchTypes_[patchi]);

            if (patchTypes_[patchi] == "symp")
            {
                patchTypes_[patchi] = "symplane";
            }
            else if (patchTypes_[patchi] == "cycl")
            {
                patchTypes_[patchi] = "cyclic";
            }
            else if (patchTypes_[patchi] == "baff")
            {
                patchTypes_[patchi] = "baffle";
            }
            else if (patchTypes_[patchi] == "moni")
            {
                patchTypes_[patchi] = "monitoring";
            }
        }

        // Create a name if needed
        if (!fndName)
        {
            patchNames_[patchi] =
                patchTypes_[patchi] + "_" + name(origRegion[patchi]);
        }
    }

    // Enforce name "Default_Boundary_Region"
    patchNames_[nPatches-1] = defaultBoundaryName;

    // Sort according to ascending region numbers, but leave
    // "Default_Boundary_Region" as the final patch
    {
        labelList sortedIndices
        (
            sortedOrder(SubList<label>(origRegion, nPatches-1))
        );

        labelList oldToNew = identity(nPatches);
        forAll(sortedIndices, i)
        {
            oldToNew[sortedIndices[i]] = i;
        }

        inplaceReorder(oldToNew, origRegion);
        inplaceReorder(oldToNew, patchTypes_);
        inplaceReorder(oldToNew, patchNames_);
        inplaceReorder(oldToNew, nPatchFaces);
    }

    // re-sort to have baffles near the end
    nBafflePatches = 1;
    if (nBafflePatches)
    {
        labelList oldToNew = identity(nPatches);
        label newIndex = 0;
        label baffleIndex = (nPatches-1 - nBafflePatches);

        for (label i=0; i < oldToNew.size()-1; ++i)
        {
            if (patchTypes_[i] == "baffle")
            {
                oldToNew[i] = baffleIndex++;
            }
            else
            {
                oldToNew[i] = newIndex++;
            }
        }

        inplaceReorder(oldToNew, origRegion);
        inplaceReorder(oldToNew, patchTypes_);
        inplaceReorder(oldToNew, patchNames_);
        inplaceReorder(oldToNew, nPatchFaces);
    }

    mapToFoamPatchId.setSize(maxId+1, -1);
    forAll(origRegion, patchi)
    {
        mapToFoamPatchId[origRegion[patchi]] = patchi;
    }

    boundaryIds_.setSize(nPatches);
    forAll(boundaryIds_, patchi)
    {
        boundaryIds_[patchi].setSize(nPatchFaces[patchi]);
        nPatchFaces[patchi] = 0;
    }


    // Pass 2:
    //
    if (nPatches > 1 && mapToFoamCellId_.size() > 1)
    {
        IFstream is(inputName);
        readHeader(is, STARCDCore::HEADER_BND);

        while (is.read(tok).good() && tok.isLabel())
        {
            // Ignore boundary id (not needed)

            is
                >> starCellId
                >> cellFaceId
                >> starRegion
                >> configNumber
                >> patchType;

            label patchi = mapToFoamPatchId[starRegion];

            // zero-based indexing
            cellFaceId--;

            label cellId = -1;

            // convert to FOAM cell number
            if (starCellId < mapToFoamCellId_.size())
            {
                cellId = mapToFoamCellId_[starCellId];
            }

            if (cellId < 0)
            {
                Info
                    << "Boundaries inconsistent with cell file. "
                    << "Star cell " << starCellId << " does not exist"
                    << endl;
            }
            else
            {
                // restrict lookup to volume cells (no baffles)
                if (cellId < cellShapes_.size())
                {
                    label mapIndex = cellShapes_[cellId].model().index();
                    if (shapeLookup.found(mapIndex))
                    {
                        mapIndex = shapeLookup[mapIndex];
                        cellFaceId =
                            STARCDCore::starToFoamFaceAddr
                            [mapIndex][cellFaceId];
                    }
                }
                else
                {
                    // we currently use cellId >= nCells to tag baffles,
                    // we can also use a negative face number
                    cellFaceId = -1;
                }

                boundaryIds_[patchi][nPatchFaces[patchi]] =
                    cellFaceIdentifier(cellId, cellFaceId);

#ifdef DEBUG_BOUNDARY
                Info<< "bnd " << cellId << " " << cellFaceId << endl;
#endif
                // Increment counter of faces in current patch
                ++nPatchFaces[patchi];
            }
        }
    }

    // retain original information in patchPhysicalTypes_ - overwrite latter
    patchPhysicalTypes_.setSize(patchTypes_.size());


    forAll(boundaryIds_, patchi)
    {
        // resize - avoid invalid boundaries
        if (nPatchFaces[patchi] < boundaryIds_[patchi].size())
        {
            boundaryIds_[patchi].setSize(nPatchFaces[patchi]);
        }

        word origType = patchTypes_[patchi];
        patchPhysicalTypes_[patchi] = origType;

        if (origType == "symplane")
        {
            patchTypes_[patchi] = symmetryPolyPatch::typeName;
            patchPhysicalTypes_[patchi] = patchTypes_[patchi];
        }
        else if (origType == "wall")
        {
            patchTypes_[patchi] = wallPolyPatch::typeName;
            patchPhysicalTypes_[patchi] = patchTypes_[patchi];
        }
        else if (origType == "cyclic")
        {
            // incorrect. should be cyclicPatch but this
            // requires info on connected faces.
            patchTypes_[patchi] = oldCyclicPolyPatch::typeName;
            patchPhysicalTypes_[patchi] = patchTypes_[patchi];
        }
        else if (origType == "baffle")
        {
            // incorrect. tag the patch until we get proper support.
            // set physical type to a canonical "baffle"
            patchTypes_[patchi] = emptyPolyPatch::typeName;
            patchPhysicalTypes_[patchi] = "baffle";
        }
        else
        {
            patchTypes_[patchi] = polyPatch::typeName;
        }

        Info<< "patch " << patchi
            << " (region " << origRegion[patchi]
            << ": " << origType << ") type: '" << patchTypes_[patchi]
            << "' name: " << patchNames_[patchi] << endl;
    }

    // cleanup
    mapToFoamCellId_.clear();
    cellShapes_.clear();
}


//
// remove unused points
//
void Foam::fileFormats::STARCDMeshReader::cullPoints()
{
    label nPoints = points_.size();
    labelList oldToNew(nPoints, -1);

    // loop through cell faces and note which points are being used
    forAll(cellFaces_, celli)
    {
        const faceList& faces = cellFaces_[celli];
        forAll(faces, i)
        {
            const labelList& labels = faces[i];
            forAll(labels, j)
            {
                oldToNew[labels[j]]++;
            }
        }
    }

    // The new ordering and the count of unused points
    label pointi = 0;
    forAll(oldToNew, i)
    {
        if (oldToNew[i] >= 0)
        {
            oldToNew[i] = pointi++;
        }
    }

    // Report unused points
    if (nPoints > pointi)
    {
        Info<< "Unused    points  = " << (nPoints - pointi) << endl;
        nPoints = pointi;

        // Adjust points and truncate
        inplaceReorder(oldToNew, points_);
        points_.setSize(nPoints);

        // Adjust cellFaces - with mesh shapes this might be faster
        for (faceList& faces : cellFaces_)
        {
            for (face& f : faces)
            {
                inplaceRenumber(oldToNew, f);
            }
        }

        // Adjust baffles
        for (face& f : baffleFaces_)
        {
            inplaceRenumber(oldToNew, f);
        }
    }
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

bool Foam::fileFormats::STARCDMeshReader::readGeometry(const scalar scaleFactor)
{
    readPoints
    (
        starFileName(geometryFile_, STARCDCore::VRT_FILE),
        scaleFactor
    );
    readCells
    (
        starFileName(geometryFile_, STARCDCore::CEL_FILE)
    );
    cullPoints();
    readBoundary
    (
        starFileName(geometryFile_, STARCDCore::BND_FILE)
    );

    return true;
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::fileFormats::STARCDMeshReader::STARCDMeshReader
(
    const fileName& prefix,
    const objectRegistry& registry,
    const scalar scaleFactor,
    const bool keepSolids
)
:
    meshReader(prefix, scaleFactor),
    keepSolids_(keepSolids),
    cellShapes_(0),
    mapToFoamPointId_(0),
    mapToFoamCellId_(0)
{
    readAux(registry);
}


// ************************************************************************* //