voxelMeshSearch.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) 2017-2023 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 "voxelMeshSearch.H"
#include "polyMesh.H"
#include "IOobject.H"
#include "fvMesh.H"
#include "block.H"
#include "emptyPolyPatch.H"
#include "processorPolyPatch.H"
#include "fvMeshTools.H"

/* * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * */

namespace Foam
{
    defineTypeNameAndDebug(voxelMeshSearch, 0);
}


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

Foam::labelVector Foam::voxelMeshSearch::offset
(
    const labelVector& nDivs
)
{
    return labelVector(1, nDivs.x(), nDivs.x()*nDivs.y());
}


Foam::label Foam::voxelMeshSearch::index
(
    const labelVector& nDivs,
    const labelVector& voxel
)
{
    return voxel.x()+voxel.y()*nDivs.x()+voxel.z()*nDivs.x()*nDivs.y();
}


Foam::labelVector Foam::voxelMeshSearch::index3
(
    const labelVector& nDivs,
    label voxeli
)
{
    const label nxy = nDivs.x()*nDivs.y();

    labelVector voxel;
    voxel.z() = voxeli/nxy;
    voxeli = voxeli % nxy;
    voxel.y() = voxeli/nDivs.x();
    voxel.x() = voxeli%nDivs.x();

    return voxel;
}


Foam::labelVector Foam::voxelMeshSearch::index3
(
    const boundBox& bb,
    const labelVector& g,
    const point& pt
)
{
    const vector s(cmptDivide(bb.span(), vector(g.x(), g.y(), g.z())));

    labelVector v
    (
        floor((pt.x()-bb.min().x())/s.x()),
        floor((pt.y()-bb.min().y())/s.y()),
        floor((pt.z()-bb.min().z())/s.z())
    );

    return v;
}


Foam::label Foam::voxelMeshSearch::index
(
    const boundBox& bb,
    const labelVector& g,
    const point& pt,
    const bool clip
)
{
    const vector s(cmptDivide(bb.span(), vector(g.x(), g.y(), g.z())));

    labelVector v
    (
        floor((pt.x()-bb.min().x())/s.x()),
        floor((pt.y()-bb.min().y())/s.y()),
        floor((pt.z()-bb.min().z())/s.z())
    );

    if (clip)
    {
        v[0] = max(0, min(g[0]-1, v[0]));
        v[1] = max(0, min(g[1]-1, v[1]));
        v[2] = max(0, min(g[2]-1, v[2]));
    }
    else if
    (
        v[0] < 0
     || v[1] < 0
     || v[2] < 0
     || v[0] >= g[0]
     || v[1] >= g[1]
     || v[2] >= g[2]
    )
    {
        return -1;
    }

    return index(g, v);
}


Foam::point Foam::voxelMeshSearch::centre
(
    const boundBox& bb,
    const labelVector& g,
    const labelVector& voxel
)
{
    const vector s(cmptDivide(bb.span(), vector(g.x(), g.y(), g.z())));

    return bb.min()+0.5*s+point(voxel[0]*s[0], voxel[1]*s[1], voxel[2]*s[2]);
}


void Foam::voxelMeshSearch::writeGrid
(
    OBJstream& os,
    const boundBox& bb,
    const labelVector& g
)
{
    const vector s(cmptDivide(bb.span(), vector(g.x(), g.y(), g.z())));

    for (label i = 1; i < g[0]; i++)
    {
        for (label j = 0; j < g[1]; j++)
        {
            for (label k = 0; k < g[2]; k++)
            {
                point p1(bb.min()+point((i-1)*s[0], j*s[1], k*s[2]));
                point p2(bb.min()+point(i*s[0], j*s[1], k*s[2]));
                os.writeLine(p1, p2);
            }
        }
    }
    for (label i = 0; i < g[0]; i++)
    {
        for (label j = 1; j < g[1]; j++)
        {
            for (label k = 0; k < g[2]; k++)
            {
                point p1(bb.min()+point(i*s[0], (j-1)*s[1], k*s[2]));
                point p2(bb.min()+point(i*s[0], j*s[1], k*s[2]));
                os.writeLine(p1, p2);
            }
        }
    }
    for (label i = 0; i < g[0]; i++)
    {
        for (label j = 0; j < g[1]; j++)
        {
            for (label k = 1; k < g[2]; k++)
            {
                point p1(bb.min()+point(i*s[0], j*s[1], (k-1)*s[2]));
                point p2(bb.min()+point(i*s[0], j*s[1], k*s[2]));
                os.writeLine(p1, p2);
            }
        }
    }
}


Foam::label Foam::voxelMeshSearch::searchProcPatch
(
    const label faceID,
    const point& searchPoint
) const
{
    const pointField& cellCentres = mesh_.cellCentres();
    const polyBoundaryMesh& bMeshes = mesh_.boundaryMesh();

    const label patchi = bMeshes.patchID(faceID);
    const polyPatch& bMeshPatch = bMeshes[patchi];

    if (!isA<processorPolyPatch>(bMeshPatch))
    {
        return -1;
    }
    else
    {
        // Find nearest cell. Linear search since cheaper than constructing
        // tree?
        const labelUList& faceCells = bMeshPatch.faceCells();
        scalar minProximity = GREAT;

        label nearestCellI = -1;
        forAll(faceCells, i)
        {
            const point& cc = cellCentres[faceCells[i]];
            scalar proximity = magSqr(cc-searchPoint);
            if (proximity < minProximity)
            {
                minProximity = proximity;
                nearestCellI = faceCells[i];
            }
        }
        return nearestCellI;
    }
}


Foam::label Foam::voxelMeshSearch::findIntersectedFace
(
    const label cellI,
    const point& p
) const
{
    // Return -1 or the label of the face intersected when tracking from
    // p to centre of cellI

    const faceList& faces = mesh_.faces();
    const pointField& faceCentres = mesh_.faceCentres();
    const pointField& points = mesh_.points();

    const point& cc = mesh_.cellCentres()[cellI];
    const labelList& cFaces = mesh_.cells()[cellI];

    const vector q(cc-p);

    forAll(cFaces, cFacei)
    {
        label facei = cFaces[cFacei];

        pointHit hitInfo = faces[facei].intersection
        (
            p,
            q,
            faceCentres[facei],
            points,
            intersection::HALF_RAY
        );

        if (hitInfo.hit() && (hitInfo.distance() < 1))
        {
            return facei;
        }
    }
    return -1;
}


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

Foam::voxelMeshSearch::voxelMeshSearch
(
    const polyMesh& mesh,
    const bool doUpdate
)
:
    mesh_(mesh)
{
    // Determine number of voxels from number of cells in mesh
    const labelVector& dim = mesh_.geometricD();

    // Guarantee at least one voxel
    label nCells = max(1, mesh_.nCells());

    label nDivs = -1;
    if (mesh_.nGeometricD() == 1)
    {
        nDivs = nCells;
    }
    else if (mesh_.nGeometricD() == 2)
    {
        nDivs = label(Foam::sqrt(scalar(nCells)));
    }
    else
    {
        nDivs = label(Foam::cbrt(scalar(nCells)));
    }

    nDivs_ = labelVector(nDivs, nDivs, nDivs);
    forAll(dim, i)
    {
        if (dim[i] == -1)
        {
            nDivs_[i] = 1;
        }
    }

    // Redo the local bounding box
    localBb_ = boundBox(mesh_.points(), false);

    const point eps(1e-10, 1e-10, 1e-10);

    localBb_.min() = localBb_.min()-eps;
    localBb_.max() = localBb_.max()+eps;

    if (debug)
    {
        Pout<< "voxelMeshSearch : mesh:" << mesh_.name()
            << " nDivs:" << nDivs_ << endl;
    }

    if (doUpdate)
    {
        update();
    }
}


Foam::voxelMeshSearch::voxelMeshSearch
(
    const polyMesh& mesh,
    const boundBox& localBb,
    const labelVector& nDivs,
    const bool doUpdate
)
:
    mesh_(mesh),
    localBb_(localBb),
    nDivs_(nDivs)
{
    if (doUpdate)
    {
        update();
    }
}


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

bool Foam::voxelMeshSearch::update()
{
    // Initialise seed cell array

    seedCell_.setSize(cmptProduct(nDivs_));
    seedCell_ = -1;


    // Find seed cells
    const pointField& points = mesh_.points();
    const labelListList& cellPoints = mesh_.cellPoints();

    forAll(cellPoints, celli)
    {
        const labelList& cPoints = cellPoints[celli];

        // Get cell bounding box
        boundBox bb(points, cPoints, false);

        fill(seedCell_, localBb_, nDivs_, bb, celli);
    }


    if (debug)
    {
        Pout<< "voxelMeshSearch : mesh:" << mesh_.name()
            << " nDivs:" << nDivs_
            << " localBb:" << localBb_ << endl;
    }


    //const pointField& cellCentres = mesh_.cellCentres();
    //forAll(cellCentres, celli)
    //{
    //    label voxeli = index(cellCentres[celli]);
    //    seedCell_[voxeli] = celli;
    //}

    return true;
}


Foam::label Foam::voxelMeshSearch::findCell(const point& p) const
{
    // First check if the point is contained in the bounding box, else exit
    if (!localBb_.contains(p))
    {
        return -1;
    }


    // Locate the voxel index for this point. Do not clip.
    label voxeli = index(localBb_, nDivs_, p, false);

    // The point may still be inside the bb but outside the actual domain.
    if (voxeli < 0)
    {
        return -1;
    }
    else
    {
        // Inverse map to compute the seed cell.
        label celli = seedCell_[voxeli];

        if (celli < 0)
        {
            return -1;
        }
        else
        {
            // Simplified, non-parallel tracking from cell centre of
            // celli to wanted location p. Note that the cell thus
            // found does not have to be the absolute 'correct' one as
            // long as at least one of the processors finds a cell.

            track_.clear();
            while (true)
            {
                if (track_.size() < 5)
                {
                    track_.append(celli);
                }

                // I am in celli now. How many faces do I have ?
                label facei = findIntersectedFace(celli, p);

                if (facei == -1)
                {
                    return celli;
                }

                const label startOfTrack(max(0, track_.size()-5));

                label nextCell;
                if (mesh_.isInternalFace(facei))
                {
                    label own = mesh_.faceOwner()[facei];
                    label nei = mesh_.faceNeighbour()[facei];
                    nextCell = (own == celli ? nei : own);

                    if (track_.found(nextCell, startOfTrack))
                    {
                        return celli;
                    }
                }
                else
                {
                    nextCell = searchProcPatch(facei, p);

                    if (nextCell == -1 || nextCell == celli)
                    {
                        return nextCell;
                    }
                    else if (track_.found(nextCell, startOfTrack))
                    {
                        return -1;  // point is really out
                    }
                }

                celli = nextCell;
            }
            return -1;
        }
    }
}


Foam::autoPtr<Foam::fvMesh> Foam::voxelMeshSearch::makeMesh
(
    const IOobject& io
) const
{
    cellShapeList cellShapes;
    faceListList boundary;
    pointField points;
    {
        //Info<< "Creating block" << endl;

        block b
        (
            cellShape(cellModel::HEX, identity(8)),
            localBb_.points(),
            blockEdgeList(),
            blockFaceList(),
            nDivs_,
            UList<gradingDescriptors>::null()
        );

        cellShapes = b.shapes();

        //Info<< "Creating boundary faces" << endl;

        boundary.setSize(b.boundaryPatches().size());
        forAll(boundary, patchi)
        {
            faceList faces(b.boundaryPatches()[patchi].size());
            forAll(faces, facei)
            {
                faces[facei] = face(b.boundaryPatches()[patchi][facei]);
            }
            boundary[patchi].transfer(faces);
        }

        points.transfer(const_cast<pointField&>(b.points()));
    }

    //Info<< "Creating patch dictionaries" << endl;
    wordList patchNames(boundary.size());
    forAll(patchNames, patchi)
    {
        patchNames[patchi] = polyPatch::defaultName(patchi);
    }

    PtrList<dictionary> boundaryDicts(boundary.size());
    forAll(boundaryDicts, patchi)
    {
        boundaryDicts.set(patchi, new dictionary());
        dictionary& patchDict = boundaryDicts[patchi];
        patchDict.add("type", emptyPolyPatch::typeName);
    }

    //Info<< "Creating polyMesh" << endl;
    IOobject polyIO(io);
    polyIO.readOpt(IOobject::NO_READ);
    polyMesh mesh
    (
        //IOobject
        //(
        //    polyMesh::defaultRegion,
        //    runTime.constant(),
        //    runTime,
        //    IOobject::NO_READ
        //),
        polyIO,
        std::move(points),
        cellShapes,
        boundary,
        patchNames,
        boundaryDicts,
        "defaultFaces",
        emptyPolyPatch::typeName,
        false
    );

    //Info<< "Writing polyMesh" << endl;
    mesh.write();

    //Info<< "Reading fvMesh" << endl;

    fvMeshTools::createDummyFvMeshFiles
    (
        io.db(),
        io.name()
    );
    IOobject fvIO(io);
    fvIO.readOpt(IOobject::MUST_READ);

    return autoPtr<fvMesh>::New(fvIO);
}


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