faceHeatShading.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2015 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/>.
 
\*---------------------------------------------------------------------------*/
 
#include "faceHeatShading.H"
#include "fvMesh.H"
#include "boundaryRadiationProperties.H"
#include "OFstream.H"
#include "cyclicAMIPolyPatch.H"
#include "volFields.H"
#include "distributedTriSurfaceMesh.H"
 
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
    defineTypeNameAndDebug(faceHeatShading, 1);
}
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 
void Foam::faceHeatShading::writeRays
(
    const fileName& fName,
    const DynamicField<point>& endCf,
    const pointField& myFc
)
{
    OFstream str(fName);
    label vertI = 0;
 
    //Pout<< "Dumping rays to " << str.name() << endl;
 
    forAll(myFc, faceI)
    {
        meshTools::writeOBJ(str, myFc[faceI]);
        vertI++;
        meshTools::writeOBJ(str, endCf[faceI]);
        vertI++;
        str << "l " << vertI-1 << ' ' << vertI << nl;
    }
    //string cmd("objToVTK " + fName + " " + fName.lessExt() + ".vtk");
    //Pout<< "cmd:" << cmd << endl;
    // system(cmd);
}
 
 
Foam::triSurface Foam::faceHeatShading::triangulate
(
    const labelHashSet& includePatches,
    const List<labelHashSet>& includeAllFacesPerPatch
)
{
    const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
 
    // Storage for surfaceMesh. Size estimate.
    DynamicList<labelledTri> triangles
    (
        mesh_.nFaces() - mesh_.nInternalFaces()
    );
 
    label newPatchI = 0;
 
    forAllConstIter(labelHashSet, includePatches, iter)
    {
        const label patchI = iter.key();
        const polyPatch& patch = bMesh[patchI];
        const pointField& points = patch.points();
 
        label nTriTotal = 0;
 
        if (includeAllFacesPerPatch[patchI].size() > 0)
        {
            forAllConstIter
            (
                labelHashSet,
                includeAllFacesPerPatch[patchI],
                iter1
            )
            {
                const label patchFaceI = iter1.key();
 
                const face& f = patch[patchFaceI];
 
                faceList triFaces(f.nTriangles(points));
 
                label nTri = 0;
 
                f.triangles(points, nTri, triFaces);
 
                forAll(triFaces, triFaceI)
                {
                    const face& f = triFaces[triFaceI];
 
                    triangles.append
                    (
                        labelledTri(f[0], f[1], f[2], newPatchI)
                    );
                    nTriTotal++;
                }
            }
            newPatchI++;
        }
    }
 
    triangles.shrink();
 
    // Create globally numbered tri surface
    triSurface rawSurface(triangles, mesh_.points());
 
    // Create locally numbered tri surface
    triSurface surface
    (
        rawSurface.localFaces(),
        rawSurface.localPoints()
    );
 
    // Add patch names to surface
    surface.patches().setSize(newPatchI);
 
    newPatchI = 0;
 
    forAllConstIter(labelHashSet, includePatches, iter)
    {
        const label patchI = iter.key();
        const polyPatch& patch = bMesh[patchI];
 
        if (includeAllFacesPerPatch[patchI].size() > 0)
        {
            surface.patches()[newPatchI].name() = patch.name();
            surface.patches()[newPatchI].geometricType() = patch.type();
 
            newPatchI++;
        }
    }
 
    return surface;
}
 
 
void Foam::faceHeatShading::calculate()
{
    const radiation::boundaryRadiationProperties& boundaryRadiation =
        radiation::boundaryRadiationProperties::New(mesh_);
 
    label nFaces = 0;          //total number of direct hit faces
 
    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
    DynamicList<point> dynCf(nFaces);
    DynamicList<label> dynFacesI;
 
    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];
        const pointField& cf = pp.faceCentres();
 
        if (!pp.coupled() && !isA<cyclicAMIPolyPatch>(pp))
        {
            const tmp<scalarField> tt =
                boundaryRadiation.transmissivity(patchI);
            const scalarField& t = tt();
            const vectorField& n = pp.faceNormals();
 
            forAll(n, faceI)
            {
                const vector nf(n[faceI]);
                if (((direction_ & nf) > 0) && (t[faceI] == 0.0))
                {
                    dynFacesI.append(faceI + pp.start());
                    dynCf.append(cf[faceI]);
                    nFaces++;
                }
            }
        }
    }
 
    label numberPotentialHits = nFaces;
 
    reduce(numberPotentialHits, sumOp<label>());
 
    Info<< "Number of 'potential' direct hits : "
        << numberPotentialHits << endl;
 
    labelList hitFacesIds(nFaces);
    hitFacesIds.transfer(dynFacesI);
 
    pointField Cfs(hitFacesIds.size());
    Cfs.transfer(dynCf);
 
    // * * * * * * * * * * * * * * *
    // Create distributedTriSurfaceMesh
    Random rndGen(653213);
 
    // Determine mesh bounding boxes:
    List<treeBoundBox> meshBb
    (
        1,
        treeBoundBox
        (
            boundBox(mesh_.points(), false)
        ).extend(rndGen, 1e-3)
    );
 
    // Dummy bounds dictionary
    dictionary dict;
    dict.add("bounds", meshBb);
    dict.add
    (
        "distributionType",
        distributedTriSurfaceMesh::distributionTypeNames_
        [
            distributedTriSurfaceMesh::FROZEN
        ]
    );
    dict.add("mergeDistance", SMALL);
 
    labelHashSet includePatches;
    List<labelHashSet> includeAllFacesPerPatch(patches.size());
 
    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];
        if (!pp.coupled() && !isA<cyclicAMIPolyPatch>(pp))
        {
            includePatches.insert(patchI);
 
            const tmp<scalarField> tt =
                boundaryRadiation.transmissivity(patchI);
            const scalarField& tau = tt();
 
            forAll(pp, faceI)
            {
                if (tau[faceI] == 0.0)
                {
                    includeAllFacesPerPatch[patchI].insert
                    (
                        faceI //pp.start()
                    );
                }
            }
        }
    }
 
    labelList triSurfaceToAgglom(5*nFaces);
 
    triSurface localSurface = triangulate
    (
        includePatches,
        includeAllFacesPerPatch
    );
 
    distributedTriSurfaceMesh surfacesMesh
    (
        IOobject
        (
            "opaqueSurface.stl",
            mesh_.time().constant(),    // directory
            "triSurface",               // instance
            mesh_.time(),               // registry
            IOobject::NO_READ,
            //IOobject::NO_WRITE
            IOobject::AUTO_WRITE
        ),
        localSurface,
        dict
    );
 
    surfacesMesh.searchableSurface::write();
 
    triSurfaceToAgglom.resize(surfacesMesh.size());
 
    scalar maxBounding = 5.0*mag(mesh_.bounds().max() - mesh_.bounds().min());
 
    reduce(maxBounding, maxOp<scalar>());
 
    // Calculate index of faces which have a direct hit (local)
//  Info << "qian zai de htiface number is " << nFaces << endl;
 
    DynamicList<label> rayStartFace(nFaces + 0.01*nFaces);
//  Info<< "raystartFace.size is "  << rayStartFace.size() << endl;
 
    // Shoot Rays
    // * * * * * * * * * * * * * * * *
    {
 
        DynamicField<point> start(nFaces);
        DynamicField<point> end(start.size());
        DynamicList<label> startIndex(start.size());
//  Info<< "start.size is "  << start.size() << endl;
//  Info<< "end.size is "  << end.size() << endl;
 
//  Info<< "startIndex.size is "  << startIndex.size() << endl;
 
        label i = 0;
//  Info << "Now cfs size is " << Cfs.size() << endl;
//  Info << "Now hitFacesIds size is " << hitFacesIds.size() << endl;
//  Info << "direction_ " << direction_<<endl;
//  Info << "maxBounding " << maxBounding<<endl;
        do
        {
            for (; i < Cfs.size(); i++)
            {
                const point& fc = Cfs[i];
 
                const label myFaceId = hitFacesIds[i];
 
                const vector d(direction_*maxBounding);
                const vector smallDistance(0.1, 0.1, 0.1);
 
                //start.append(fc + SMALL*d);
                start.append(fc - 0.1*direction_);
//Info << "fc-samll*d is " << fc-SMALL*d << endl;
//Info << "fc-samll*d is " << fc-smallDistance << endl;
//Info << "myFaceId is " << myFaceId << endl;
 
 
                startIndex.append(myFaceId);
 
                end.append(fc - d);
//Info << "fc-d is " << fc-d << endl<<endl;
 
            }
 
        }while (returnReduce(i < Cfs.size(), orOp<bool>()));
//Info <<"start number is" <<start.size() << endl;
//Info <<"end number is" <<end.size() << endl;
//Info <<"startIndex size number is" <<startIndex.size() << endl;
 
        List<pointIndexHit> hitInfo(startIndex.size());
//Info <<"=====before findline HitInfo size number is " <<hitInfo.size() << endl<<endl;
forAll (hitInfo,rayI)
{
//  Info <<rayI << " now hit yes or no is " << hitInfo[rayI].hit()<< endl;
}
 
        surfacesMesh.findLine(start, end, hitInfo);
        //surfacesMesh.findLineAll(start, end, hitInfo);
//Info <<"=-=====after fineline HitInfo size number is " <<hitInfo.size() << endl<<endl;
forAll (hitInfo,rayI)
{
    //Info <<rayI << " now hit yes or no is " << hitInfo[rayI].hit()<< endl;
if (hitInfo[rayI].hit())
{
    //Info <<rayI << " now hit point is " << hitInfo[rayI].hitPoint()<< endl;
}
else
{
    //Info <<rayI << " now hit miss point is " << hitInfo[rayI].missPoint()<< endl;
}
 
    //Info <<rayI << " now hit index is " << hitInfo[rayI].index()<< endl;
//hitInfo[rayI].write(Info);
//Info <<endl<<endl;
 
}
 
        // Collect the rays which has 'only one not wall' obstacle bettween
        // start and end.
        // If the ray hit itself get stored in dRayIs
        forAll (hitInfo, rayI)
        {
            if (!hitInfo[rayI].hit())
            {
//          Info << hitInfo[rayI].hit() << endl;
//          Info << hitInfo[rayI].rawPoint() << endl;
//          hitInfo[rayI].write(Info) ;
//          Info <<rayI <<  "ray can not hit wall ,blocked by other wall" << endl;
        //  hitInfo[rayI].write(Info) ;
            //Info<<endl <<rayI <<  " ray can not hit wall ,blocked by other wall" << endl;
    //      Info<<endl << "startIndex[rayi] is "<< startIndex[rayI]<< endl;
            rayStartFace.append(startIndex[rayI]);
    
    //      Info<< endl<< rayI  <<" rayStartFace size Is "<< rayStartFace.size() << endl;
            }
            else
            {
        //  hitInfo[rayI].write(Info) ;
            //Info<<endl <<rayI <<  " ray can not hit wall ,blocked by other wall" << endl;
            //Info<< endl <<"rayStartFace size Is "<< rayStartFace.size() << endl;
            //Info << hitInfo[rayI].hit() << endl;
            //Info << hitInfo[rayI].rawPoint() << endl;
            //hitInfo[rayI].write(Info) ;
            //rayStartFace.append(startIndex[rayI]);
    //      Info<<endl <<rayI <<  " ray can not hit wall ,blocked by other wall" << endl;
    //      Info<<rayI  <<" rayStartFace size Is "<< rayStartFace.size() << endl;
            }
        //  Info << endl;
        }
 
//Info << "now rayStartFace size is " << rayStartFace.size()<<endl;
 
        // Plot all rays between visible faces.
        if (debug)
        {
            writeRays
            (
                mesh_.time().path()/"allVisibleFaces.obj",
                end,
                Cfs
            );
        }
 
        start.clear();
        startIndex.clear();
        end.clear();
    }
 
    rayStartFaces_.transfer(rayStartFace);
 
    if (debug)
    {
        tmp<volScalarField> thitFaces
        (
            new volScalarField
            (
                IOobject
                (
                    "hitFaces",
                    mesh_.time().timeName(),
                    mesh_,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                mesh_,
                dimensionedScalar("zero", dimless, 0)
            )
        );
 
        volScalarField& hitFaces = thitFaces();
Info << "now hitFace size is " << hitFaces.size()<<endl;
 
        hitFaces.boundaryField() = 0.0;
        forAll (rayStartFaces_, i)
        {
            const label faceI = rayStartFaces_[i];
            label patchID = patches.whichPatch(faceI);
            const polyPatch& pp = patches[patchID];
            hitFaces.boundaryField()[patchID][faceI - pp.start()] = 1.0;
        }
//Info << "after for loop, now hitFace size is " << hitFaces.size()<<endl;
//Info << "write hit faces" << endl;
//下面这一句执行的话，会生成1文件夹，然后里面只有一个hitface，没有其他变量。
//hitFaces.write();
    }
 
    label totalHitFaces = rayStartFaces_.size();
 
    reduce(totalHitFaces, sumOp<label>());
 
    Info<< "Total number of hit faces : " <<  totalHitFaces << endl;
}
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::faceHeatShading::faceHeatShading
(
    const fvMesh& mesh,
    const vector dir,
    const labelList& hitFaceList
)
:
    mesh_(mesh),
    direction_(dir),
    rayStartFaces_(hitFaceList)
{}
 
 
 
Foam::faceHeatShading::faceHeatShading
(
    const fvMesh& mesh,
    const vector dir
)
:
    mesh_(mesh),
    direction_(dir),
    rayStartFaces_(0)
{
    calculate();
}
 
 
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
Foam::faceHeatShading::~faceHeatShading()
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
void Foam::faceHeatShading::correct()
{
    calculate();
}
 
 
// ************************************************************************* //