ENH: Added new createViewFactors utility

Creates view factors for the view factor radiation model.

User-selectable models:

- raySearchEngine: model to generate rays, i.e. face-to-face connections
- viewFactorModel: model to compute the view factors

For visualisation, use:

- Write the view factors as a volume field

    writeViewFactors    yes;

- Write the rays using OBJ format:

    writeRays       yes; // default = no

Participating patches must be in the \c vewFactorWall group, i.e. using the
\c inGroups entry of the "\<case\>/polyMesh/boundary" file.

\verbatim
myPatch
{
    type            wall;
    inGroups        2(wall viewFactorWall);
    ...
}
\endverbatim

Reads:

- <constant>/viewFactorsDict : main controls
- <constant>/finalAgglom : agglomeration addressing (from faceAgglomerate)

Generates:

- <constant>/F : view factors (matrix)
- <constant>/mapDist : map used for parallel running
- <constant>/globalFaceFaces : face addressing
This commit is contained in:
Andrew Heather
2024-06-14 16:53:28 +01:00
committed by Kutalmis Bercin
parent 7c45670c8b
commit ab5f6dbf41
23 changed files with 4115 additions and 0 deletions

View File

@ -0,0 +1,11 @@
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments
#------------------------------------------------------------------------------
wmake $targetType viewFactorModels
wmake $targetType createViewFactors
#------------------------------------------------------------------------------

View File

@ -0,0 +1,3 @@
createViewFactors.C
EXE = $(FOAM_APPBIN)/createViewFactors

View File

@ -0,0 +1,7 @@
EXE_INC = \
-I../viewFactorModels/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lviewFactorModels \
-lfiniteVolume

View File

@ -0,0 +1,165 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 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/>.
Application
createViewFactors
Group
grpPreProcessingUtilities
Description
Creates view factors to be used in the view-factor radiation model.
Operands:
\table
Operand | Type | Location
input | dictionary | \<constant\>/viewFactorsDict
input | dictionary | \<constant\>/finalAgglom
output | scalarListList | \<constant\>/F
output | mapDistribute | \<constant\>/mapDist
output | labelListList | \<constant\>/globalFaceFaces
output | volScalarField | \<time\>/viewVectorField
output | OBJ | allVisibleFaces.obj
\endtable
where the dictionaries mean:
\table
Dictionary | Description
viewFactorsDict | Main-control dictionary
finalAgglom | (Optional) Agglomeration addressing (from faceAgglomerate)
F | View factors (matrix)
mapDist | Map used for parallel running
globalFaceFaces | Face addressing
viewVectorField | View factors as a volume field
allVisibleFaces.obj | The visualisation of the rays
\endtable
Usage
Minimal example in \c <constant>/viewFactorsDict:
\verbatim
// Inherited entries
raySearchEngine <word>;
agglomerate <bool>;
nRayPerFace <label>;
writeViewFactors <bool>;
writeRays <bool>;
...
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
raySearchEngine | Ray search engine type | word | yes | -
agglomerate | Flag to agglomeration | bool | yes | -
nRayPerFace | Number of rays issued per face | label | yes | -
writeViewFactors | Flag to write the view factor field | bool | yes |-
writeRays | Flag to write the ray geometry | bool | no | false
\endtable
Options for the \c raySearchEngine entry:
\verbatim
voxel | Ray search engine discretising space into uniform voxels
\endverbatim
The inherited entries are elaborated in:
- \link viewFactorModel.H \endlink
- \link raySearchEngine.H \endlink
Note
- Participating patches must be in the \c vewFactorWall group, i.e. using the
\c inGroups entry of the "\<case\>/polyMesh/boundary" file.
\verbatim
myPatch
{
type wall;
inGroups 2(wall viewFactorWall);
...
}
\endverbatim
Reads:
- <constant>/viewFactorsDict : main controls
- <constant>/finalAgglom : agglomeration addressing (from faceAgglomerate)
Generates:
- <constant>/F : view factors (matrix)
- <constant>/mapDist : map used for parallel running
- <constant>/globalFaceFaces : face addressing
SeeAlso
- Foam::VF::raySearchEngine
- Foam::VE::viewFactorModel
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "Time.H"
#include "fvMesh.H"
#include "viewFactorModel.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "addRegionOption.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createNamedMesh.H"
IOdictionary dict
(
IOobject
(
"viewFactorsDict",
runTime.constant(),
mesh,
IOobject::MUST_READ
)
);
// Calculate the view factors
auto modelPtr = VF::viewFactorModel::New(mesh, dict);
modelPtr->calculate();
Info<< nl;
runTime.printExecutionTime(Info);
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,12 @@
raySearchEngine/raySearchEngine/raySearchEngine.C
raySearchEngine/raySearchEngine/raySearchEngineNew.C
raySearchEngine/voxel/voxelRaySearchEngine.C
/* raySearchEngine/nonUniformVoxel/nonUniformVoxelRaySearchEngine.C */
/* raySearchEngine/allToAll/allToAllRaySearchEngine.C */
viewFactorModel/viewFactorModel/viewFactorModel.C
viewFactorModel/viewFactorModel/viewFactorModelNew.C
viewFactorModel/viewFactor2AI/viewFactor2AI.C
viewFactorModel/viewFactor2LI/viewFactor2LI.C
viewFactorModel/viewFactorHottel/viewFactorHottel.C
LIB = $(FOAM_LIBBIN)/libviewFactorModels

View File

@ -0,0 +1,14 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/parallel/distributed/lnInclude \
-I$(LIB_SRC)/fileFormats/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lsurfMesh \
-lmeshTools \
-ldistributed \
-lradiationModels \
-lfileFormats

View File

@ -0,0 +1,620 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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 "raySearchEngine.H"
#include "surfaceFields.H"
#include "volFields.H"
#include "meshTools.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace VF
{
defineTypeNameAndDebug(raySearchEngine, 0);
defineRunTimeSelectionTable(raySearchEngine, mesh);
const label raySearchEngine::maxDynListLength = 1000000000;
}
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::VF::raySearchEngine::check
(
const labelList& nVisibleFaceFaces
)
{
label nRay = 0;
label nFaceMin = labelMax;
label nFaceMax = labelMin;
for (const label n : nVisibleFaceFaces)
{
nFaceMin = min(nFaceMin, n);
nFaceMax = max(nFaceMax, n);
nRay += n;
}
const label nFace = nVisibleFaceFaces.size();
const label nGlobalRays = returnReduce(nRay, sumOp<label>());
if (nGlobalRays == 0)
{
FatalErrorInFunction
<< "No rays identified - view factors will not be calculated"
<< exit(FatalError);
}
const label globalNFacesMin = returnReduce(nFaceMin, minOp<label>());
const label globalNFacesMax = returnReduce(nFaceMax, maxOp<label>());
const label nGlobalFaces = returnReduce(nFace, sumOp<label>());
const scalar avgFaces = nGlobalRays/scalar(nGlobalFaces);
Info<< "\nRay summary:" << nl
<< " Number of rays: " << nGlobalRays << nl
<< " Number of rays-per-face (min, max, average): ("
<< globalNFacesMin << ", "
<< globalNFacesMax << ", "
<< avgFaces << ")" << endl;
}
Foam::label Foam::VF::raySearchEngine::closestPointIndex
(
const point& p0,
const List<point>& pts
)
{
label pointi = -1;
scalar distSqr = GREAT;
forAll(pts, pti)
{
const scalar d2 = magSqr(pts[pti] - p0);
if (d2 < distSqr)
{
pointi = pti;
distSqr = d2;
}
}
return pointi;
}
void Foam::VF::raySearchEngine::createAgglomeration(const IOobject& io)
{
Info<< "\nFace agglomeration: active" << nl
<< " Reading file " << io.name() << endl;
// Read agglomeration map
const labelListIOList finalAgglom(io);
Info<< " Creating coarse mesh" << nl;
agglomMeshPtr_.reset
(
new singleCellFvMesh
(
IOobject
(
IOobject::scopedName("agglom", mesh_.name()),
mesh_.time().timeName(),
mesh_.time(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
finalAgglom
)
);
const auto& coarseMesh = agglomMeshPtr_();
// Calculate total number of fine and coarse faces
nCoarseFace_ = 0;
nFace_ = 0;
const polyBoundaryMesh& finePatches = mesh_.boundaryMesh();
const polyBoundaryMesh& coarsePatches = coarseMesh.boundaryMesh();
for (const label patchi : patchIDs_)
{
nCoarseFace_ += coarsePatches[patchi].size();
nFace_ += finePatches[patchi].size();
}
Info<< "\nTotal number of coarse faces: "
<< returnReduce(nCoarseFace_, sumOp<label>())
<< endl;
Info<< "\nTotal number of fine faces: "
<< returnReduce(nFace_, sumOp<label>())
<< endl;
// Collect local Cf, Sf, agglom index on coarse mesh
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
DynamicList<point> localCf(nCoarseFace_);
DynamicList<vector> localSf(nCoarseFace_);
DynamicList<label> localAgg(nCoarseFace_);
for (const label patchi : patchIDs_)
{
const labelList& agglom = finalAgglom[patchi];
if (agglom.empty()) continue;
label nAgglom = max(agglom) + 1;
const labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
const labelList& coarsePatchFace = coarseMesh.patchFaceMap()[patchi];
const pointField& coarseCf = coarseMesh.Cf().boundaryField()[patchi];
const vectorField& coarseSf = coarseMesh.Sf().boundaryField()[patchi];
const polyPatch& pp = finePatches[patchi];
patchAreas_[patchi] += sum(coarseMesh.magSf().boundaryField()[patchi]);
forAll(coarseCf, facei)
{
const label coarseFacei = coarsePatchFace[facei];
const label agglomi = coarseFacei + coarsePatches[patchi].start();
// Construct single coarse face
const labelList& fineFaces = coarseToFine[coarseFacei];
uindirectPrimitivePatch cf
(
UIndirectList<face>(pp, fineFaces),
pp.points()
);
// Collect all points (vertices, face centres)
const label nFaces = cf.faceCentres().size();
const label nPoints = cf.localPoints().size();
List<point> allPoints(nFaces + nPoints);
SubList<point>(allPoints, nFaces) = cf.faceCentres();
SubList<point>(allPoints, nPoints, nFaces) = cf.localPoints();
// Coarse face centre set to closest point
const label pti = closestPointIndex(coarseCf[facei], allPoints);
if (pti != -1)
{
localCf.push_back(allPoints[pti]);
localSf.push_back(coarseSf[facei]);
localAgg.push_back(agglomi);
}
}
}
Info<< "\nAssembled coarse patch data" << endl;
// Distribute local coarse Cf and Sf for shooting rays
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
allCf_[Pstream::myProcNo()].transfer(localCf);
allSf_[Pstream::myProcNo()].transfer(localSf);
allAgg_[Pstream::myProcNo()].transfer(localAgg);
Pstream::allGatherList(allCf_);
Pstream::allGatherList(allSf_);
Pstream::allGatherList(allAgg_);
Pstream::listCombineGather(patchAreas_, plusEqOp<scalar>());
Pstream::broadcast(patchAreas_);
globalNumbering_ = globalIndex(nCoarseFace_);
}
void Foam::VF::raySearchEngine::createGeometry()
{
DynamicList<point> localCf(mesh_.nBoundaryFaces());
DynamicList<vector> localSf(mesh_.nBoundaryFaces());
const auto& pbm = mesh_.boundaryMesh();
for (const label patchi : patchIDs_)
{
localCf.push_back(pbm[patchi].faceCentres());
localSf.push_back(pbm[patchi].faceAreas());
patchAreas_[patchi] += sum(mesh_.magSf().boundaryField()[patchi]);
}
Info<< "\nAssembled patch data" << endl;
nFace_ = localCf.size();
nCoarseFace_ = -1;
allCf_[Pstream::myProcNo()].transfer(localCf);
allSf_[Pstream::myProcNo()].transfer(localSf);
Pstream::allGatherList(allCf_);
Pstream::allGatherList(allSf_);
Pstream::listCombineGather(patchAreas_, plusEqOp<scalar>());
Pstream::broadcast(patchAreas_);
globalNumbering_ = globalIndex(nFace_);
}
void Foam::VF::raySearchEngine::createParallelAddressing
(
labelList& rayEndFace
) const
{
// Construct compact numbering
// - return map from remote to compact indices
// (per processor (!= myProcNo) a map from remote index to compact index)
// - construct distribute map
// - renumber rayEndFace into compact addressing
DebugInfo << "\nCreating map distribute" << endl;
List<Map<label>> compactMap(Pstream::nProcs());
mapPtr_.reset(new mapDistribute(globalNumbering_, rayEndFace, compactMap));
DebugInfo << "\nCreating compact-to-global addressing" << endl;
// Invert compactMap (from processor+localface to compact) to go
// from compact to processor+localface (expressed as a globalIndex)
compactToGlobal_.resize_nocopy(mapPtr_->constructSize());
// Local indices first
// Note: are not in compactMap
for (label i = 0; i < globalNumbering_.localSize(); ++i)
{
compactToGlobal_[i] = globalNumbering_.toGlobal(i);
}
forAll(compactMap, proci)
{
const Map<label>& localToCompactMap = compactMap[proci];
forAllConstIters(localToCompactMap, iter)
{
compactToGlobal_[*iter] =
globalNumbering_.toGlobal(proci, iter.key());
}
}
}
Foam::coordSystem::cartesian Foam::VF::raySearchEngine::createCoordSystem
(
const point& origin,
const vector& dir
) const
{
vector axis(Zero);
for (direction d=0; d<3; ++d)
{
axis = dir^tensor::I.col(d);
// Remove empty direction for 2D
if (mesh_.nSolutionD() == 2)
{
meshTools::constrainDirection(mesh_, mesh_.solutionD(), axis);
}
if (magSqr(axis) > 0)
{
axis.normalise();
break;
}
}
return coordSystem::cartesian(origin, dir, axis);
}
Foam::tmp<Foam::pointField> Foam::VF::raySearchEngine::createHemiPoints
(
const label nRayPerFace
) const
{
auto themiPts = tmp<pointField>::New(nRayPerFace);
auto& hemiPts = themiPts.ref();
const label nPoints = hemiPts.size();
if (mesh_.nSolutionD() == 3)
{
// Point in range -1 < x < 1; -1 < y < 1; 0 < z 1
forAll(hemiPts, pointi)
{
const scalar phi = Foam::acos(1 - (pointi + 0.5)/nPoints);
const scalar theta =
mathematical::pi*(1 + Foam::sqrt(5.0))*(pointi + 0.5);
hemiPts[pointi] =
vector
(
Foam::cos(theta)*Foam::sin(phi),
Foam::sin(theta)*Foam::sin(phi),
Foam::cos(phi)
);
}
}
else if (mesh_.nSolutionD() == 2)
{
// Point in range -1 < x < 1; y = 0; 0 < z < 1; _\|/_
forAll(hemiPts, pointi)
{
const scalar theta = mathematical::pi*(pointi+0.5)/nPoints;
hemiPts[pointi] = vector(Foam::cos(theta), 0, Foam::sin(theta));
}
}
return themiPts;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::VF::raySearchEngine::raySearchEngine
(
const fvMesh& mesh,
const dictionary& dict
)
:
mesh_(mesh),
mapPtr_(nullptr),
compactToGlobal_(),
globalNumbering_(),
patchGroup_(dict.getOrDefault<word>("patchGroup", "viewFactorWall")),
patchIDs_(mesh_.boundaryMesh().indices(patchGroup_)),
patchAreas_(mesh_.boundaryMesh().nNonProcessor(), Zero),
agglomerate_(dict.get<bool>("agglomerate")),
agglomMeshPtr_(nullptr),
nFace_(-1),
nCoarseFace_(-1),
allCf_(UPstream::nProcs()),
allSf_(UPstream::nProcs()),
allAgg_(UPstream::nProcs())
{
Info<< "\nParticipating patches:" << endl;
forAll(patchIDs_, i)
{
const label patchi = patchIDs_[i];
Info<< " " << i << ": " << mesh_.boundaryMesh()[patchi].name()
<< endl;
}
const word agglomName(dict.getOrDefault<word>("agglom", "finalAgglom"));
IOobject agglomIO
(
agglomName,
mesh_.facesInstance(),
mesh_,
IOobject::MUST_READ
);
if (agglomerate_)
{
// Sets allCf_, allSf_, allAgg_ based on coarse mesh
// Sets nFace_, nCoarseFace_
createAgglomeration(agglomIO);
}
else
{
// Check for presence of finalAgglom - will cause problems in later
// calculations with viewFactor radiation model
if (agglomIO.typeHeaderOk<labelListIOList>())
{
WarningInFunction
<< "Found agglomeration file: " << agglomIO.objectPath() << nl
<< " This is inconsistent with the view factor calculation "
<< "and should be removed" << nl << endl;
}
// Sets allCf_, allSf_ based on fine mesh
// Sets nFace_; nCoarseFace_ remains unset (-1)
createGeometry();
}
globalNumbering_ =
nCoarseFace_ == -1 ? globalIndex(nFace_) : globalIndex(nCoarseFace_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::VF::raySearchEngine::correct
(
labelListList& visibleFaceFaces
) const
{
labelList rayStartFace;
labelList rayEndFace;
shootRays(rayStartFace, rayEndFace);
const label nFace = nParticipatingFaces();
// Calculate number of visible faces from each local start face
labelList nVisibleFaceFaces(nFace, Zero);
for (const label facei : rayStartFace)
{
++nVisibleFaceFaces[facei];
}
check(nVisibleFaceFaces);
createParallelAddressing(rayEndFace);
// Set visible face-faces
// visibleFaceFaces has:
// (local face, local viewed face) = compact viewed face
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
visibleFaceFaces.resize_nocopy(nFace);
forAll(nVisibleFaceFaces, facei)
{
visibleFaceFaces[facei].resize_nocopy(nVisibleFaceFaces[facei]);
}
nVisibleFaceFaces = 0;
forAll(rayStartFace, i)
{
const label facei = rayStartFace[i];
const label sloti = rayEndFace[i];
visibleFaceFaces[facei][nVisibleFaceFaces[facei]++] = sloti;
}
}
void Foam::VF::raySearchEngine::compactAddressing
(
const mapDistribute& map,
pointField& compactCf,
vectorField& compactSf,
List<List<vector>>& compactFineSf,
List<List<point>>& compactFineCf,
DynamicList<List<point>>& compactPoints,
DynamicList<label>& compactPatchId
) const
{
compactCf.resize_nocopy(map.constructSize());
compactSf.resize_nocopy(map.constructSize());
compactFineSf.resize_nocopy(map.constructSize());
compactFineCf.resize_nocopy(map.constructSize());
compactPoints.setCapacity(map.constructSize());
compactPatchId.setCapacity(map.constructSize());
// Insert my local values area and centre values
if (agglomMeshPtr_)
{
SubList<vector>(compactSf, nCoarseFace_) = allSf_[Pstream::myProcNo()];
SubList<point>(compactCf, nCoarseFace_) = allCf_[Pstream::myProcNo()];
const auto& coarseMesh = agglomMeshPtr_();
const auto& coarsePatches = coarseMesh.boundaryMesh();
const auto& coarseFaces = coarseMesh.faces();
const auto& coarsePoints = coarseMesh.points();
const auto& finalAgglom = coarseMesh.patchFaceAgglomeration();
// Insert my fine local values per coarse face
label sloti = 0;
for (const label patchi : patchIDs_)
{
const auto& fineCfp = mesh_.Cf().boundaryField()[patchi];
const auto& fineSfp = mesh_.Sf().boundaryField()[patchi];
const labelList& agglom = finalAgglom[patchi];
if (agglom.empty()) continue;
const label nAgglom = max(agglom) + 1;
const labelListList coarseToFine = invertOneToMany(nAgglom, agglom);
const labelList& coarsePatchFace =
coarseMesh.patchFaceMap()[patchi];
const label coarseStart = coarsePatches[patchi].start();
forAll(coarseToFine, coarsei)
{
compactPatchId.push_back(patchi);
const vectorField localPoints
(
coarsePoints,
coarseFaces[coarseStart + coarsei]
);
compactPoints.push_back(localPoints);
const label coarseFacei = coarsePatchFace[coarsei];
const labelList& fineFaces = coarseToFine[coarseFacei];
List<point>& fineCf = compactFineCf[sloti];
fineCf.resize_nocopy(fineFaces.size());
fineCf = UIndirectList<point>(fineCfp, fineFaces);
List<vector>& fineSf = compactFineSf[sloti];
fineSf.resize_nocopy(fineFaces.size());
fineSf = UIndirectList<vector>(fineSfp, fineFaces);
++sloti;
}
}
}
else
{
SubList<vector>(compactSf, nFace_) = allSf_[Pstream::myProcNo()];
SubList<point>(compactCf, nFace_) = allCf_[Pstream::myProcNo()];
const auto& patches = mesh_.boundaryMesh();
const faceList& faces = mesh_.faces();
label sloti = 0;
for (const label patchi : patchIDs_)
{
const auto& Sfp = mesh_.Sf().boundaryField()[patchi];
const auto& Cfp = mesh_.Cf().boundaryField()[patchi];
const polyPatch& pp = patches[patchi];
forAll(pp, facei)
{
compactPatchId.push_back(patchi);
const auto& fpts = faces[facei + pp.start()];
compactPoints.push_back(List<point>(mesh_.points(), fpts));
compactFineCf[sloti] = List<point>({Cfp[facei]});
compactFineSf[sloti] = List<vector>({Sfp[facei]});
++sloti;
}
}
}
// Do all swapping
map.distribute(compactSf);
map.distribute(compactCf);
map.distribute(compactFineCf);
map.distribute(compactFineSf);
map.distribute(compactPoints);
map.distribute(compactPatchId);
}
// ************************************************************************* //

View File

@ -0,0 +1,289 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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/>.
Class
Foam::VF::raySearchEngine
Description
Base class for ray search engines
Participating patches must be in the \c viewFactorWall group, i.e. using the
\c inGroups entry of the "\<case\>/polyMesh/boundary" file.
\verbatim
myPatch
{
type wall;
inGroups 2(wall viewFactorWall);
...
}
\endverbatim
Face agglomeration can be employed, created using the \c faceAgglomerate
utility. The file name to be read can be user-defined:
\verbatim
// Name of agglomeration file; default = finalAgglom
agglom finalAgglom;
\endverbatim
SourceFiles
raySearchEngine.C
\*---------------------------------------------------------------------------*/
#ifndef Foam_vf_raySearchEngine_H
#define Foam_vf_raySearchEngine_H
#include "cartesianCS.H"
#include "mapDistribute.H"
#include "singleCellFvMesh.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace VF
{
/*---------------------------------------------------------------------------*\
Class raySearchEngine Declaration
\*---------------------------------------------------------------------------*/
class raySearchEngine
{
protected:
// Protected Data
//- Reference to the mesh
const fvMesh& mesh_;
//- Parallel map
mutable autoPtr<mapDistribute> mapPtr_;
//- Compact to global addressing
mutable labelList compactToGlobal_;
//- Global numbering
globalIndex globalNumbering_;
//- Name of patch group to identify participating patches
const word patchGroup_;
//- List of participating patch IDs
labelList patchIDs_;
//- Patch areas
scalarList patchAreas_;
//- Agglomeration flag
bool agglomerate_;
//- Agglomerated mesh representation
autoPtr<singleCellFvMesh> agglomMeshPtr_;
//- Number of original faces
label nFace_;
//- Number of coarse faces
label nCoarseFace_;
//- List of all face centres per processor
List<pointField> allCf_;
//- List of all face areas per processor
List<vectorField> allSf_;
//- List of all face agglomeration index per processor
List<labelField> allAgg_;
// Protected Member Functions
static void check(const labelList& nVisibleFaceFaces);
static label closestPointIndex
(
const point& p0,
const List<point>& pts
);
//- Create patch geometry based on the original mesh
void createGeometry();
//- Create parallel addressing - map, compact-to-global
void createParallelAddressing(labelList& rayEndFace) const;
//- Create Cartesian co-ordinate system
coordSystem::cartesian createCoordSystem
(
const point& origin,
const vector& dir
) const;
//- Create patch geometry based on the agglomerated mesh
void createAgglomeration(const IOobject& io);
//- Create a set of points describing a hemisphere
// Note: origin is (0 0 0)
tmp<pointField> createHemiPoints(const label nRayPerFace) const;
public:
static const label maxDynListLength;
//- Run-time type information
TypeName("raySearchEngine");
//- Selection table
declareRunTimeSelectionTable
(
autoPtr,
raySearchEngine,
mesh,
(
const fvMesh& mesh,
const dictionary& dict
),
(mesh, dict)
);
//- Selector
static autoPtr<raySearchEngine> New
(
const fvMesh& mesh,
const dictionary& dict
);
// Generated Methods
//- No copy construct
raySearchEngine(const raySearchEngine&) = delete;
//- No copy assignment
void operator=(const raySearchEngine&) = delete;
//- Constructor
raySearchEngine(const fvMesh& mesh, const dictionary& dict);
//- Destructor
virtual ~raySearchEngine() = default;
// Public Member Functions
// Access
//- Reference to the mesh
inline const fvMesh& mesh() const noexcept;
//- Parallel map
inline const mapDistribute& map() const;
//- Compact to global addressing
inline const labelList& compactToGlobal() const noexcept;
//- Global numbering
inline const globalIndex& globalNumbering() const noexcept;
//- List of participating patch IDs
inline const labelList& patchIDs() const noexcept;
//- Patch areas
inline const scalarList& patchAreas() const noexcept;
//- Number of participating faces
inline label nParticipatingFaces() const;
//- List of all face centres per processor
inline const List<pointField>& allCf() const noexcept;
//- List of all face areas per processor
inline const List<vectorField>& allSf() const noexcept;
//- List of all face agglomeration index per processor
inline const List<labelField>& allAgg() const noexcept;
// Main calculation functions
//- Shoot rays; returns lists of ray start and end faces
virtual void shootRays
(
labelList& rayStartFaceOut,
labelList& rayEndFaceOut
) const = 0;
//- Correct
virtual void correct(labelListList& visibleFaceFaces) const;
//- Create compact addressing
void compactAddressing
(
const mapDistribute& map,
pointField& compactCf,
vectorField& compactSf,
List<List<vector>>& compactFineSf,
List<List<point>>& compactFineCf,
DynamicList<List<point>>& compactPoints,
DynamicList<label>& compactPatchId
) const;
//- Interpolate field
template<class Type>
void interpolate
(
GeometricField<Type, fvPatchField, volMesh>& fld,
const List<List<Type>>& values
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace VF
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "raySearchEngineI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "raySearchEngineTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,101 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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/>.
\*---------------------------------------------------------------------------*/
const Foam::fvMesh& Foam::VF::raySearchEngine::mesh() const noexcept
{
return mesh_;
}
const Foam::mapDistribute& Foam::VF::raySearchEngine::map() const
{
if (!mapPtr_)
{
FatalErrorInFunction
<< "mapDistribute has not been set"
<< abort(FatalError);
}
return mapPtr_;
}
const Foam::labelList&
Foam::VF::raySearchEngine::compactToGlobal() const noexcept
{
return compactToGlobal_;
}
const Foam::globalIndex&
Foam::VF::raySearchEngine::globalNumbering() const noexcept
{
return globalNumbering_;
}
const Foam::labelList& Foam::VF::raySearchEngine::patchIDs() const noexcept
{
return patchIDs_;
}
const Foam::scalarList& Foam::VF::raySearchEngine::patchAreas() const noexcept
{
return patchAreas_;
}
Foam::label Foam::VF::raySearchEngine::nParticipatingFaces() const
{
if (nCoarseFace_ == -1) return nFace_;
return nCoarseFace_;
}
const Foam::List<Foam::pointField>&
Foam::VF::raySearchEngine::allCf() const noexcept
{
return allCf_;
}
const Foam::List<Foam::vectorField>&
Foam::VF::raySearchEngine::allSf() const noexcept
{
return allSf_;
}
const Foam::List<Foam::labelField>&
Foam::VF::raySearchEngine::allAgg() const noexcept
{
return allAgg_;
}
// ************************************************************************* //

View File

@ -0,0 +1,59 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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 "raySearchEngine.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::VF::raySearchEngine> Foam::VF::raySearchEngine::New
(
const fvMesh& mesh,
const dictionary& dict
)
{
const word modelType(dict.get<word>("raySearchEngine"));
Info<< "Selecting " << typeName << ": " << modelType << endl;
auto* ctorPtr = meshConstructorTable(modelType);
if (!ctorPtr)
{
FatalIOErrorInLookup
(
dict,
typeName,
modelType,
*meshConstructorTablePtr_
) << exit(FatalIOError);
}
return autoPtr<raySearchEngine>(ctorPtr(mesh, dict));
}
// ************************************************************************* //

View File

@ -0,0 +1,86 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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/>.
\*---------------------------------------------------------------------------*/
template<class Type>
void Foam::VF::raySearchEngine::interpolate
(
GeometricField<Type, fvPatchField, volMesh>& fld,
const List<List<Type>>& values
) const
{
label compacti = 0;
auto& vfbf = fld.boundaryFieldRef();
if (agglomMeshPtr_)
{
const auto& coarseMesh = agglomMeshPtr_();
const auto& finalAgglom = coarseMesh.patchFaceAgglomeration();
for (const label patchi : patchIDs_)
{
const labelList& agglom = finalAgglom[patchi];
if (agglom.empty()) continue;
label nAgglom = max(agglom) + 1;
const labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
const labelList& coarsePatchFace =
coarseMesh.patchFaceMap()[patchi];
forAll(coarseToFine, i)
{
const label coarseFacei = coarsePatchFace[i];
const labelList& fineFaces = coarseToFine[coarseFacei];
const Type sumValues = sum(values[compacti]);
for (const label fineFacei : fineFaces)
{
vfbf[patchi][fineFacei] = sumValues;
}
++compacti;
}
}
}
else
{
label compacti = 0;
for (const label patchi : patchIDs_)
{
auto& vfp = vfbf[patchi];
for (Type& vfi : vfp)
{
vfi = sum(values[compacti++]);
}
}
}
}
// ************************************************************************* //

View File

@ -0,0 +1,948 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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 "voxelRaySearchEngine.H"
#include "processorPolyPatch.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace VF
{
defineTypeNameAndDebug(voxel, 0);
addToRunTimeSelectionTable(raySearchEngine, voxel, mesh);
}
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::VF::voxel::setTriangulation(const fvMesh& mesh)
{
Info<< "\nCreating triangulated surface" << endl;
// Storage for surfaceMesh. Size estimate.
DynamicList<labelledTri> triangles(mesh.nBoundaryFaces());
DynamicList<label> globalFaces(mesh.nBoundaryFaces());
label nFace = 0;
const auto& pbm = mesh.boundaryMesh();
forAll(patchIDs_, i)
{
const label patchi = patchIDs_[i];
const polyPatch& patch = pbm[patchi];
const pointField& points = patch.points();
for (const face& f : patch)
{
label nTri = 0;
faceList triFaces(f.nTriangles(points));
f.triangles(points, nTri, triFaces);
const label globalFacei =
globalNumbering_.toGlobal(Pstream::myProcNo(), nFace++);
for (const face& f : triFaces)
{
triangles.push_back(labelledTri(f[0], f[1], f[2], i));
globalFaces.push_back(globalFacei);
}
}
}
triToGlobalFace_.transfer(globalFaces);
Info<< " Total number of triangles: "
<< returnReduce(triangles.size(), sumOp<label>())
<< endl;
triangles.shrink();
surface_ = triSurface(triangles, mesh.points());
surface_.compactPoints();
}
Foam::labelList Foam::VF::voxel::setFaceVertexHits
(
const fvMesh& mesh,
const labelList& patchIDs
)
{
labelList vertHits(mesh.points().size(), Zero);
if (mesh.nSolutionD() == 3)
{
const auto& pbm = mesh.boundaryMesh();
// Create a new triangulation based on the surface agglomeration
for (const label patchI : patchIDs)
{
const polyPatch& patch = pbm[patchI];
for (const face& f : patch)
{
for (const label fpi : f)
{
++vertHits[fpi];
}
}
}
for (const auto& pp : pbm)
{
const labelList& meshPts = pp.meshPoints();
if (pp.size())
{
if (isA<processorPolyPatch>(pp))
{
// Add all processor patch points
for (const label pi : meshPts)
{
++vertHits[pi];
}
}
else
{
// Add boundary points
const auto& bndyPts = pp.boundaryPoints();
for (const label pi : bndyPts)
{
++vertHits[meshPts[pi]];
}
}
}
}
}
return vertHits;
}
void Foam::VF::voxel::setCoarseTriangulation(const fvMesh& mesh)
{
Info<< "\nCreating triangulated surface" << endl;
// Filter out fine mesh points along coarse mesh faces
// Storage for surfaceMesh. Size estimate.
DynamicList<labelledTri> triangles(mesh.nBoundaryFaces());
DynamicList<label> globalFaces(mesh.nBoundaryFaces());
labelList vertHits = setFaceVertexHits(mesh, patchIDs_);
// Only simplifying edges for 3D
const label nVertMin = mesh.nSolutionD() == 3 ? 2 : 0;
label nInvalid = 0;
label nFace = 0;
const auto& pbm = mesh.boundaryMesh();
for (const label patchi : patchIDs_)
{
const polyPatch& patch = pbm[patchi];
const pointField& points = patch.points();
for (const face& f : patch)
{
DynamicList<label> faceVerts;
for (const label fpi : f)
{
if (vertHits[fpi] > nVertMin)
{
faceVerts.push_back(fpi);
}
}
if (faceVerts.size() < 3)
{
++nInvalid;
continue;
}
label nTri = 0;
const face reducedFace(faceVerts);
faceList triFaces(reducedFace.nTriangles(points));
reducedFace.triangles(points, nTri, triFaces);
const label globalFacei =
globalNumbering_.toGlobal(Pstream::myProcNo(), nFace++);
for (const face& f : triFaces)
{
triangles.push_back(labelledTri(f[0], f[1], f[2], patchi));
globalFaces.push_back(globalFacei);
}
}
}
triToGlobalFace_.transfer(globalFaces);
Info<< " Total number of triangles: "
<< returnReduce(triangles.size(), sumOp<label>())
<< "\n Number of invalid (removed) triangles: "
<< returnReduce(nInvalid, sumOp<label>())
<< endl;
triangles.shrink();
surface_ = triSurface(triangles, mesh.points());
surface_.compactPoints();
}
void Foam::VF::voxel::broadcast()
{
// Every processor has whole surface
const globalIndex globalFaceIdx(globalIndex::gatherOnly{}, surface_.size());
const globalIndex globalPointIdx
(
globalIndex::gatherOnly{},
surface_.points().size()
);
List<labelledTri> globalSurfaceTris(globalFaceIdx.gather(surface_));
pointField globalSurfacePoints(globalPointIdx.gather(surface_.points()));
List<label> globalTriToGlobalFace(globalFaceIdx.gather(triToGlobalFace_));
for (const label proci : globalPointIdx.allProcs())
{
const label offset = globalPointIdx.localStart(proci);
if (offset)
{
for
(
labelledTri& tri
: globalSurfaceTris.slice(globalFaceIdx.range(proci))
)
{
tri[0] += offset;
tri[1] += offset;
tri[2] += offset;
}
}
}
surface_ =
triSurface
(
std::move(globalSurfaceTris),
std::move(globalSurfacePoints)
);
Pstream::broadcast(surface_);
triToGlobalFace_ = std::move(globalTriToGlobalFace);
Pstream::broadcast(triToGlobalFace_);
}
Foam::pointHit Foam::VF::voxel::rayTriIntersect
(
const label trii,
const point& origin,
const vector& direction
) const
{
// Note: origin was made relative to voxel mesh on entry to hit function
// - need to convert back into global position to be consistent with
// triangles for intersection tests
const auto& ind = surface_[trii];
const auto& pts = surface_.points();
// Note: flipped orientation of triangle (into domain) so that we can use
// visibility check when performing ray-triangle intersections
const triPointRef tri(pts[ind[0]], pts[ind[2]], pts[ind[1]]);
static scalar tolerance = 1e-6;
return
tri.intersection
(
globalPosition(origin),
direction,
intersection::VISIBLE,
tolerance
);
}
void Foam::VF::voxel::writeBox
(
OBJstream& os,
bool lines,
const boundBox& bb
) const
{
os.write(treeBoundBox(bb), lines);
}
void Foam::VF::voxel::writeVoxels(const word& fName) const
{
if (!UPstream::master()) return;
OBJstream os(fName);
Info<< "Writing voxels to " << os.name() << endl;
boundBox bb;
const bool lines = true;
for (label i = 0; i < nijk_[0]; ++i)
{
for (label j = 0; j < nijk_[1]; ++j)
{
for (label k = 0; k < nijk_[2]; ++k)
{
bb.min() = voxelMin(i, j, k);
bb.max() = voxelMax(i, j, k);
writeBox(os, lines, bb);
}
}
}
Info<< "- done" << endl;
}
void Foam::VF::voxel::writeTriBoundBoxes(const word& fName) const
{
if (!UPstream::master()) return;
OBJstream os(fName);
Info<< "Writing triangle boundBoxes to " << os.name() << endl;
const bool lines = true;
for (const auto& voxeli : objects_)
{
for (const label trii : voxeli)
{
writeBox(os, lines, objectBbs_[trii]);
}
}
Info<< "- done" << endl;
}
void Foam::VF::voxel::writeHitObjects
(
const label voxeli,
const point& origin,
const vector& dir
) const
{
OBJstream os("voxel_" + Foam::name(voxeli) + ".obj");
// Write voxel
labelVector ijk = this->ijk(voxeli);
boundBox voxelBb;
voxelBb.min() = voxelMin(ijk[0], ijk[1], ijk[2]);
voxelBb.max() = voxelMax(ijk[0], ijk[1], ijk[2]);
writeBox(os, true, voxelBb);
for (const auto& trii : objects_[voxeli])
{
writeBox(os, true, objectBbs_[trii]);
const auto& ind = surface_[trii];
const auto& pts = surface_.points();
const triPointRef tri(pts[ind[0]], pts[ind[1]], pts[ind[2]]);
os.write(tri, false);
}
}
Foam::pointIndexHit Foam::VF::voxel::hitObject
(
const label voxeli,
const point& origin,
const vector& dir,
const vector& t,
const scalar minDistance
) const
{
if (objects_[voxeli].empty()) return pointIndexHit();
// boundBox rayBb;
// rayBb.add(origin);
// rayBb.add(origin + dir*(dir & t));
label triHiti = -1;
//rayBb.add(origin + dir);
//rayBb.inflate(0.01);
if (debug > 2)
{
writeHitObjects(voxeli, origin, dir);
}
// Determine all triangles that intersect with ray
// - only keep nearest hit
pointHit closestHit;
for (const auto& trii : objects_[voxeli])
{
// Only perform ray/tri intersection if bound boxes intersect
//if (objectBbs_[trii].overlaps(rayBb))
{
pointHit pi = rayTriIntersect(trii, origin, dir);
if
(
pi.hit()
&& (
pi.distance() > minDistance
&& pi.distance() < closestHit.distance()
)
)
{
triHiti = trii;
closestHit = pi;
}
}
}
return pointIndexHit(closestHit, triHiti);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::VF::voxel::voxel(const fvMesh& mesh, const dictionary& dict)
:
raySearchEngine(mesh, dict),
bb0_(),
span0_(Zero),
nijk_(Zero),
dxyz_(Zero),
nRayPerFace_(dict.get<label>("nRayPerFace")),
nTriPerVoxelMax_(dict.getOrDefault<label>("nTriPerVoxelMax", 50)),
depthMax_(dict.getOrDefault<label>("depthMax", 5)),
objects_(),
objectBbs_()
{
if (agglomMeshPtr_)
{
setCoarseTriangulation(*agglomMeshPtr_);
}
else
{
setTriangulation(mesh);
}
broadcast();
objectBbs_.resize_nocopy(surface_.size());
bb0_.add(surface_.points());
bb0_.inflate(0.01);
span0_ = bb0_.span();
{
scalar maxDx = span0_.x();
scalar maxDy = span0_.y();
scalar maxDz = span0_.z();
const label refDn = 8;
scalar maxDim = max(maxDx, max(maxDy, maxDz));
setVoxelDims
(
refDn*maxDx/maxDim,
refDn*maxDy/maxDim,
refDn*maxDz/maxDim
);
objects_.resize_nocopy(nVoxel());
}
label depth = 0;
label trii = 0;
voxelise(objects_, trii, depth);
Info<< "\nCreated voxel mesh: " << nijk_ << endl;
if ((debug > 3) && UPstream::master())
{
writeVoxels("voxels.obj");
writeTriBoundBoxes("triBoundBoxes.obj");
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::VF::voxel::refineObjects
(
List<DynamicList<label>>& objects,
const label triMax
)
{
refineVoxelDims();
if (debug > 2) Pout<< "Refining voxels: n=" << nijk_ << endl;
List<DynamicList<label>> objectsNew(objects.size()*8);
for (label trii = 0; trii <= triMax; ++trii)
{
addBbToVoxels(objectBbs_[trii], trii, objectsNew);
}
objects.transfer(objectsNew);
}
Foam::label Foam::VF::voxel::addBbToVoxels
(
const boundBox& bb,
const label trii,
List<DynamicList<label>>& objects
) const
{
//Pout<< "addBbToVoxels: trii=" << trii << " n=" << nijk_ << endl;
const point minbb(localPosition(bb.min()));
const label i0 = max(0, floor(minbb.x()/dxyz_[0]));
const label j0 = max(0, floor(minbb.y()/dxyz_[1]));
const label k0 = max(0, floor(minbb.z()/dxyz_[2]));
const point maxbb(localPosition(bb.max()));
const label i1 = min(nijk_[0], ceil(maxbb.x()/dxyz_[0]));
const label j1 = min(nijk_[1], ceil(maxbb.y()/dxyz_[1]));
const label k1 = min(nijk_[2], ceil(maxbb.z()/dxyz_[2]));
label nTriMax = 0;
for (label ii = i0; ii < i1; ++ii)
{
for (label jj = j0; jj < j1; ++jj)
{
for (label kk = k0; kk < k1; ++kk)
{
const label voxeli = this->voxeli(ii, jj, kk);
objects[voxeli].push_back(trii);
nTriMax = max(nTriMax, objects[voxeli].size());
}
}
}
return nTriMax;
}
void Foam::VF::voxel::voxelise
(
List<DynamicList<label>>& objects,
const label trii0,
const label depth
)
{
if (debug > 2)
{
Pout<< "voxelise - start at tri=" << trii0
<< " depth=" << depth
<< endl;
}
const auto& points = surface_.points();
for (label trii = trii0; trii < surface_.size(); ++trii)
{
// Triangle bounding box
boundBox bb(points, surface_[trii]);
bb.inflate(0.01);
objectBbs_[trii] = bb;
const label nVoxelMax = addBbToVoxels(bb, trii, objects);
// Number of triangles per voxel - if exceed limit, refine voxels...
if (nVoxelMax > nTriPerVoxelMax_ && depth < depthMax_)
{
refineObjects(objects, trii);
voxelise(objects, trii + 1, depth + 1);
break;
}
}
}
Foam::pointIndexHit Foam::VF::voxel::hit
(
const label tri0,
const vector& direction0
) const
{
if (tri0 > surface_.size()-1)
{
FatalErrorInFunction
<< "Index tri0 out of bounds: " << tri0
<< ". Surface size: " << surface_.size()
<< abort(FatalError);
}
return hit(surface_.faceCentres()[tri0], direction0);
}
Foam::pointIndexHit Foam::VF::voxel::hit
(
const point& origin0,
const vector& direction0
) const
{
// Initialise return value
pointIndexHit hitInfo;
const point origin(localPosition(origin0));
if (cmptMin(origin) < 0)
{
FatalErrorInFunction
<< "Point located outside voxel mesh"
<< " - possible coarsening problem?"
<< abort(FatalError);
}
if (magSqr(direction0) < ROOTVSMALL)
{
WarningInFunction
<< "Supplied direction has zero size"
<< endl;
return hitInfo;
}
const vector dir(normalised(direction0));
labelVector ijk(Zero);
labelVector step(Zero);
vector tDelta(vector::max);
vector tMax(vector::max);
for (direction d = 0; d < 3; ++d)
{
ijk[d] = floor(origin[d]/dxyz_[d]);
step[d] = sign0(dir[d]);
if (step[d] != 0)
{
tDelta[d] = mag(dxyz_[d]/dir[d]);
scalar voxelMax = (1 + ijk[d] - neg(dir[d]))*dxyz_[d];
tMax[d] = (voxelMax - origin[d])/dir[d];
}
}
if (debug > 2)
{
Info<< "surfBb:" << boundBox(surface_.points())
<< " bb:" << bb0_
<< " origin" << origin0
<< " voxel_origin:" << origin
<< " ijk:" << ijk
<< " step:" << step
<< " dxyz:" << dxyz_
<< " tDelta:" << tDelta
<< " tMax:" << tMax
<< endl;
}
auto traverse = [&](const label i)
{
ijk[i] += step[i];
if (outOfBounds(ijk, i)) return false;
tMax[i] += tDelta[i];
return true;
};
while (true)
{
const label voxeli = this->voxeli(ijk);
if (debug > 2)
{
Info<< "ijk:" << ijk
<< " voxeli:" << voxeli
<< " t:" << tMax
<< " objs:" << objects_[voxeli].size()
<< endl;
}
hitInfo = hitObject(voxeli, origin, dir, tMax);
if (hitInfo.hit())
{
// Valid hit
break;
}
else
{
if (tMax[0] < tMax[1] && tMax[0] < tMax[2])
{
if (!traverse(0)) break;
}
else if (tMax[1] < tMax[2])
{
if (!traverse(1)) break;
}
else
{
if (!traverse(2)) break;
}
}
}
return hitInfo;
}
void Foam::VF::voxel::shootRays
(
labelList& rayStartFaceOut,
labelList& rayEndFaceOut
) const
{
(void)mesh_.time().cpuTimeIncrement();
const pointField& myFc = allCf_[Pstream::myProcNo()];
const vectorField& myArea = allSf_[Pstream::myProcNo()];
const label nTotalRays = myFc.size()*nRayPerFace_;
if (nTotalRays > maxDynListLength)
{
FatalErrorInFunction
<< "Dynamic list needs more capacity to support " << nRayPerFace_
<< " rays per face (" << nTotalRays << "). "
<< "Limit set by maxDynListLength: " << maxDynListLength
<< abort(FatalError);
}
DynamicList<label> rayStartFace(nTotalRays);
DynamicList<label> rayEndFace(rayStartFace.size());
DynamicList<label> endFacei(nTotalRays);
DynamicList<label> startFacei(nTotalRays);
const pointField hemi(createHemiPoints(nRayPerFace_));
Info<< "\nShooting rays" << endl;
const scalar nDiv = 10;
const label nStep = floor(myFc.size()/nDiv);
labelHashSet localFaceHits;
for (label stepi=0; stepi<nDiv; ++stepi)
{
const label step0 = stepi*nStep;
const label step1 = stepi == nDiv - 1 ? myFc.size() : (stepi+1)*nStep;
for (label i = step0; i < step1; ++i)
{
// Info<< "i/N = " << i << "/" << (myFc.size()-1)
// << " step0:" << step0 << " step1:" << step1
// << endl;
localFaceHits.clear();
const point& origin = myFc[i];
const vector dir(-normalised(myArea[i]));
const coordSystem::cartesian cs = createCoordSystem(origin, dir);
const vectorField pts(cs.transformPoint(hemi));
for (const auto& p : pts)
{
const pointIndexHit hitInfo = hit(origin, p-origin);
if (hitInfo.hit())
{
label hitFacei = triToGlobalFace_[hitInfo.index()];
if (localFaceHits.insert(hitFacei))
{
endFacei.push_back(hitFacei);
startFacei.push_back(i);
}
}
else
{
// Miss
}
}
}
const label globalStepi = returnReduce(stepi, minOp<label>()) + 1;
Info<< " " << globalStepi/nDiv*100 << "% complete" << endl;
}
// Ensure all rays are unique/filter out duplicates
// - add symmetric connections for non-self-intersecting rays
if (UPstream::parRun())
{
edgeHashSet uniqueRays;
List<DynamicList<label>> remoteStartFace(Pstream::nProcs());
List<DynamicList<label>> remoteEndFace(Pstream::nProcs());
const labelList globalStartFaceIDs
(
globalNumbering_.toGlobal(startFacei)
);
forAll(globalStartFaceIDs, rayi)
{
label start = globalStartFaceIDs[rayi];
label end = endFacei[rayi];
const label endProci = globalNumbering_.whichProcID(end);
if (start > end) Swap(start, end);
if (uniqueRays.insert(edge(start, end)))
{
if (endProci != UPstream::myProcNo())
{
// Convert local face into global face and vice-versa
remoteStartFace[endProci].push_back(start);
remoteEndFace[endProci].push_back(end);
}
}
}
PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
// Send remote data
for (const int domain : Pstream::allProcs())
{
if (domain != Pstream::myProcNo())
{
UOPstream str(domain, pBufs);
str << remoteStartFace[domain]
<< remoteEndFace[domain];
}
}
pBufs.finishedSends();
// Consume
for (const int domain : Pstream::allProcs())
{
if (domain != Pstream::myProcNo())
{
UIPstream is(domain, pBufs);
is >> remoteStartFace[domain]
>> remoteEndFace[domain];
forAll(remoteStartFace[domain], i)
{
const label start = remoteStartFace[domain][i];
const label end = remoteEndFace[domain][i];
uniqueRays.insert(edge(start, end));
}
}
}
// Populate ray addressing based on uniqueRays
for (const edge& e : uniqueRays)
{
const label start = e.first();
const label end = e.second();
bool sameFace = start == end;
if (globalNumbering_.isLocal(start))
{
// Ray originates from this processor
const label localStart = globalNumbering_.toLocal(start);
rayStartFace.append(localStart);
rayEndFace.append(end);
}
if (!sameFace && globalNumbering_.isLocal(end))
{
// Ray hits this processor
// - add symmetric ray from end->start
const label localEnd = globalNumbering_.toLocal(end);
rayStartFace.append(localEnd);
rayEndFace.append(start);
}
}
}
else
{
// Populate ray addressing based on uniqueRays
edgeHashSet uniqueRays;
forAll(startFacei, rayi)
{
label start = startFacei[rayi];
label end = endFacei[rayi];
if (start > end) Swap(start, end);
if (uniqueRays.insert(edge(start, end)))
{
rayStartFace.append(start);
rayEndFace.append(end);
if (start != end)
{
// Add symmetric ray from end->start
rayStartFace.append(end);
rayEndFace.append(start);
}
}
}
}
rayStartFaceOut.transfer(rayStartFace);
rayEndFaceOut.transfer(rayEndFace);
Info<< " Time taken: " << mesh_.time().cpuTimeIncrement() << " s"
<< endl;
}
// ************************************************************************* //

View File

@ -0,0 +1,284 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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/>.
Class
Foam::VF::voxel
Description
Ray search engine based on discretising space into uniform voxels
Voxels are refined using 2x2x2 refinement.
The number of rays per face is supplied by the user, whereby rays are
issued uniformly across a hemisphere.
Usage
Minimal example by using \c <constant>/viewFactorsDict:
\verbatim
// Mandatory entries
nRayPerFace <label>;
// Optional entries
nTriPerVoxelMax <label>;
depthMax <label>;
// Inherited entries
...
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
nRayPerFace | Number of rays issued per face | label | yes | -
nRayPerFace | Target number of triangles per voxel | label | no | 50
depthMax | Maximum voxel refinement depth | label | no | 5
\endtable
The inherited entries are elaborated in:
- \link raySearchEngine.H \endlink
SourceFiles
voxelRaySearchEngine.C
SeeAlso
- Foam::VF::raySearchEngine
\*---------------------------------------------------------------------------*/
#ifndef Foam_vf_voxelRaySearchEngine_H
#define Foam_vf_voxelRaySearchEngine_H
#include "raySearchEngine.H"
#include "labelVector.H"
#include "OBJstream.H"
#include "pointIndexHit.H"
#include "triSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace VF
{
/*---------------------------------------------------------------------------*\
Class voxel Declaration
\*---------------------------------------------------------------------------*/
class voxel
:
public raySearchEngine
{
// Private Data
//- Triangulation of view factor patches
triSurface surface_;
//- Triangle to mesh face addressing
labelList triToGlobalFace_;
//- Surface bounding box
boundBox bb0_;
//- Span of surface bounding box
vector span0_;
//- Voxel discretisation
labelVector nijk_;
//- Voxel dimensions
vector dxyz_;
//- Number of rays issued per face
const label nRayPerFace_;
//- Maximum number of triangles per voxel (trigger to refine voxels)
const label nTriPerVoxelMax_;
//- Maximum depth of voxel refinement. Note: voxels are refined 2x2x2
const label depthMax_;
//- List of triangle per voxel
List<DynamicList<label>> objects_;
//- List of triangle bounding boxes
List<boundBox> objectBbs_;
// Private Member Functions
inline bool outOfBounds(const labelVector& ijk, const label dir) const;
inline point localPosition(const vector& globalPosition) const;
inline point globalPosition(const vector& localPosition) const;
inline void setVoxelDims(const label i, const label j, const label k);
inline void refineVoxelDims();
inline point voxelMin
(
const label i,
const label j,
const label k
) const;
inline point voxelMax
(
const label i,
const label j,
const label k
) const;
inline constexpr label sign0(const scalar x) const;
//- Set triangulation based on original mesh
void setTriangulation(const fvMesh& mesh);
//- Set the participating face vertices when simplifying edges
static labelList setFaceVertexHits
(
const fvMesh& mesh,
const labelList& patchIDs
);
//- Set triangulation based on agglomeration
void setCoarseTriangulation(const fvMesh& mesh);
//- Broadcast triangulation across all procs
void broadcast();
void refineObjects
(
List<DynamicList<label>>& objects,
const label triMax
);
label addBbToVoxels
(
const boundBox& bb,
const label trii,
List<DynamicList<label>>& objects
) const;
void voxelise
(
List<DynamicList<label>>& objects,
const label trii0,
const label depth
);
pointHit rayTriIntersect
(
const label trii,
const point& origin,
const vector& direction
) const;
pointIndexHit hitObject
(
const label voxeli,
const point& origin,
const vector& direction,
const vector& t,
const scalar minDistance = 1e-6
) const;
void writeBox
(
OBJstream& os,
bool lines,
const boundBox& bb
) const;
void writeVoxels(const word& fName) const;
void writeTriBoundBoxes(const word& fName) const;
void writeHitObjects
(
const label voxeli,
const point& origin,
const vector& dir
) const;
public:
//- Runtime type information
TypeName("voxel");
//- Constructor
voxel(const fvMesh& mesh, const dictionary& dict);
//- Destructor
virtual ~voxel() = default;
// Public Member Functions
inline labelVector nijk() const noexcept;
inline label nVoxel() const noexcept;
inline label voxeli(const labelVector ijk) const noexcept;
inline label voxeli
(
const label i,
const label j,
const label k
) const noexcept;
inline labelVector ijk(const label voxeli) const noexcept;
pointIndexHit hit(const point& origin, const vector& dir) const;
pointIndexHit hit(const label tri0, const vector& dir) const;
virtual void shootRays
(
labelList& rayStartFaceOut,
labelList& rayEndFaceOut
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace VF
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "voxelRaySearchEngineI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,155 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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/>.
\*---------------------------------------------------------------------------*/
bool Foam::VF::voxel::outOfBounds
(
const labelVector& ijk,
const label dir
) const
{
return (ijk[dir] < 0 || ijk[dir] >= nijk_[dir]);
};
Foam::point Foam::VF::voxel::localPosition(const vector& globalPosition) const
{
return globalPosition - bb0_.min();
}
Foam::point Foam::VF::voxel::globalPosition(const vector& localPosition) const
{
return bb0_.min() + localPosition;
}
void Foam::VF::voxel::setVoxelDims(const label i, const label j, const label k)
{
nijk_[0] = max(1, i);
nijk_[1] = max(1, j);
nijk_[2] = max(1, k);
dxyz_[0] = span0_[0]/nijk_[0];
dxyz_[1] = span0_[1]/nijk_[1];
dxyz_[2] = span0_[2]/nijk_[2];
}
void Foam::VF::voxel::refineVoxelDims()
{
nijk_ *= 2;
// Do not refine empty direction for 2D
const auto& solutionD = mesh_.solutionD();
for (direction d=0; d<3; ++d)
{
if (solutionD[d] == -1)
{
nijk_[d] = 1;
}
}
setVoxelDims(nijk_[0], nijk_[1], nijk_[2]);
}
Foam::point Foam::VF::voxel::voxelMin
(
const label i,
const label j,
const label k
) const
{
return point(i*dxyz_[0], j*dxyz_[1], k*dxyz_[2]);
}
Foam::point Foam::VF::voxel::voxelMax
(
const label i,
const label j,
const label k
) const
{
return voxelMin(i+1, j+1, k+1);
}
constexpr Foam::label Foam::VF::voxel::sign0(const scalar x) const
{
if (x > 0) return 1;
if (x < 0) return -1;
return 0;
};
Foam::labelVector Foam::VF::voxel::nijk() const noexcept
{
return nijk_;
}
Foam::label Foam::VF::voxel::nVoxel() const noexcept
{
return nijk_[0]*nijk_[1]*nijk_[2];
}
Foam::label Foam::VF::voxel::voxeli
(
const labelVector ijk
) const noexcept
{
return voxeli(ijk[0], ijk[1], ijk[2]);
}
Foam::label Foam::VF::voxel::voxeli
(
const label i,
const label j,
const label k
) const noexcept
{
return i + (nijk_[0]*(j + (nijk_[1]*k)));
}
Foam::labelVector Foam::VF::voxel::ijk(const label voxeli) const noexcept
{
const label nx = nijk_[0];
const label ny = nijk_[1];
const label i = voxeli%nx;
const label k = voxeli/nx/ny;
const label j = (voxeli/nx)%ny;
return labelVector(i, j, k);
}
// ************************************************************************* //

View File

@ -0,0 +1,138 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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 "viewFactor2AI.H"
#include "mathematicalConstants.H"
#include "addToRunTimeSelectionTable.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace VF
{
defineTypeNameAndDebug(viewFactor2AI, 0);
addToRunTimeSelectionTable(viewFactorModel, viewFactor2AI, mesh);
}
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
Foam::scalar Foam::VF::viewFactor2AI::calculateFij
(
const point& xi,
const point& xj,
const vector& dAi,
const vector& dAj
)
{
const vector r(xj - xi);
const scalar rMag = mag(r);
const scalar dAiMag(mag(dAi));
const scalar dAjMag(mag(dAj));
if (rMag > SMALL && dAiMag > SMALL && dAjMag > SMALL)
{
const vector nr(r/rMag);
const vector ni(dAi/dAiMag);
const vector nj(dAj/dAjMag);
const scalar Fij =
-(nr & ni)*(nr & nj)/sqr(rMag)*dAjMag*dAiMag/mathematical::pi;
if (Fij > 0) return Fij;
}
return 0;
}
Foam::scalarListList Foam::VF::viewFactor2AI::calculate
(
const labelListList& visibleFaceFaces,
const pointField& compactCf,
const vectorField& compactSf,
const List<List<vector>>& compactFineSf,
const List<List<point>>& compactFineCf,
const DynamicList<List<point>>& compactPoints,
const DynamicList<label>& compactPatchId
) const
{
// Fill local view factor matrix
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalarListList Fij(visibleFaceFaces.size());
forAll(visibleFaceFaces, facei)
{
if (debug > 1)
{
Pout<< "facei:" << facei << "/" << visibleFaceFaces.size()
<< endl;
}
const labelList& visibleFaces = visibleFaceFaces[facei];
Fij[facei].resize(visibleFaces.size(), Zero);
const point& dCi = compactCf[facei];
const vector& Ai = compactSf[facei];
const scalar magAi = mag(Ai);
if (magAi < SMALL) continue;
forAll(visibleFaces, visibleFacei)
{
const label sloti = visibleFaces[visibleFacei];
const point& dCj = compactCf[sloti];
const vector& Aj = compactSf[sloti];
const scalar Fij2AI = calculateFij(dCi, dCj, Ai, Aj);
Fij[facei][visibleFacei] = Fij2AI/magAi;
}
}
return Fij;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::VF::viewFactor2AI::viewFactor2AI
(
const fvMesh& mesh,
const dictionary& dict
)
:
viewFactorModel(mesh, dict)
{}
// ************************************************************************* //

View File

@ -0,0 +1,117 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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/>.
Class
Foam::VF::viewFactor2AI
Description
Computes view factors according to the double area integral (2AI) method.
Usage
Minimal example in \c <constant>/viewFactorsDict:
\verbatim
// Inherited entries
...
\endverbatim
The inherited entries are elaborated in:
- \link viewFactorModel.H \endlink
SourceFiles
viewFactorModel.C
\*---------------------------------------------------------------------------*/
#ifndef Foam_vf_viewFactor2AI_H
#define Foam_vf_viewFactor2AI_H
#include "viewFactorModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace VF
{
/*---------------------------------------------------------------------------*\
Class viewFactor2AI Declaration
\*---------------------------------------------------------------------------*/
class viewFactor2AI
:
public viewFactorModel
{
protected:
// Protected Member Functions
//- Calculate view factor using the double-area integral
static scalar calculateFij
(
const point& xi,
const point& xj,
const vector& dAi,
const vector& dAj
);
//- Calculate
virtual scalarListList calculate
(
const labelListList& visibleFaceFaces,
const pointField& compactCf,
const vectorField& compactSf,
const List<List<vector>>& compactFineSf,
const List<List<point>>& compactFineCf,
const DynamicList<List<point>>& compactPoints,
const DynamicList<label>& compactPatchId
) const;
public:
//- Runtime type information
TypeName("viewFactor2AI");
//- Constructor
viewFactor2AI(const fvMesh& mesh, const dictionary& dict);
//- Destructor
virtual ~viewFactor2AI() = default;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace VF
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,141 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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 "viewFactor2LI.H"
#include "mathematicalConstants.H"
#include "addToRunTimeSelectionTable.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace VF
{
defineTypeNameAndDebug(viewFactor2LI, 0);
addToRunTimeSelectionTable(viewFactorModel, viewFactor2LI, mesh);
}
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
Foam::scalar Foam::VF::viewFactor2LI::calculateFij
(
const List<point>& lPoints,
const List<point>& rPoints,
const scalar alpha
)
{
scalar Fij = 0;
forAll(lPoints, i)
{
// Edge vector and centroid of edge i
const vector si(lPoints[i] - lPoints.rcValue(i));
const point ci(0.5*(lPoints[i] + lPoints.rcValue(i)));
forAll(rPoints, j)
{
// Edge vector and centroid of edge j
const vector sj(rPoints[j] - rPoints.rcValue(j));
const point cj(0.5*(rPoints[j] + rPoints.rcValue(j)));
vector r(ci - cj);
if (mag(r) < SMALL)
{
r = alpha*si;
}
Fij += (si & sj)*Foam::log(r & r);
}
}
return max(0, 0.25*Fij/mathematical::pi);
}
Foam::scalarListList Foam::VF::viewFactor2LI::calculate
(
const labelListList& visibleFaceFaces,
const pointField& compactCf,
const vectorField& compactSf,
const List<List<vector>>& compactFineSf,
const List<List<point>>& compactFineCf,
const DynamicList<List<point>>& compactPoints,
const DynamicList<label>& compactPatchId
) const
{
// Fill local view factor matrix
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalarListList Fij(visibleFaceFaces.size());
forAll(visibleFaceFaces, facei)
{
if (debug > 1)
{
Pout<< "facei:" << facei << "/" << visibleFaceFaces.size()
<< endl;
}
const labelList& visibleFaces = visibleFaceFaces[facei];
Fij[facei].resize(visibleFaces.size(), Zero);
const vector& Ai = compactSf[facei];
const scalar magAi = mag(Ai);
forAll(visibleFaces, visibleFacei)
{
const label sloti = visibleFaces[visibleFacei];
const List<point>& lPoints = compactPoints[facei];
const List<point>& rPoints = compactPoints[sloti];
const scalar Fij2LI = calculateFij(lPoints, rPoints, alpha_);
Fij[facei][visibleFacei] = Fij2LI/magAi;
}
}
return Fij;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::VF::viewFactor2LI::viewFactor2LI
(
const fvMesh& mesh,
const dictionary& dict
)
:
viewFactorModel(mesh, dict),
alpha_(dict.getOrDefault("alpha", 0.21))
{}
// ************************************************************************* //

View File

@ -0,0 +1,130 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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/>.
Class
Foam::VF::viewFactor2LI
Description
Computes view factors according to the double line integral (2LI) method.
Usage
Minimal example in \c <constant>/viewFactorsDict:
\verbatim
// Optional entries
alpha <scalar>;
// Inherited entries
...
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
alpha | Perturbation for common edges | scalar | no | 0.21
\endtable
The inherited entries are elaborated in:
- \link viewFactorModel.H \endlink
SourceFiles
viewFactorModel.C
\*---------------------------------------------------------------------------*/
#ifndef Foam_vf_viewFactor2LI_H
#define Foam_vf_viewFactor2LI_H
#include "viewFactorModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace VF
{
/*---------------------------------------------------------------------------*\
Class viewFactor2LI Declaration
\*---------------------------------------------------------------------------*/
class viewFactor2LI
:
public viewFactorModel
{
protected:
// Protected Data
//- Perturbation for common edges; default = 0.21
const scalar alpha_;
// Protected Member Functions
//- Calculate view factor using the double-area integral
static scalar calculateFij
(
const List<point>& lPoints,
const List<point>& rPoints,
const scalar alpha
);
//- Calculate
virtual scalarListList calculate
(
const labelListList& visibleFaceFaces,
const pointField& compactCf,
const vectorField& compactSf,
const List<List<vector>>& compactFineSf,
const List<List<point>>& compactFineCf,
const DynamicList<List<point>>& compactPoints,
const DynamicList<label>& compactPatchId
) const;
public:
//- Runtime type information
TypeName("viewFactor2LI");
//- Constructor
viewFactor2LI(const fvMesh& mesh, const dictionary& dict);
//- Destructor
virtual ~viewFactor2LI() = default;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace VF
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,153 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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 "viewFactorHottel.H"
#include "mathematicalConstants.H"
#include "fvMesh.H"
#include "meshTools.H"
//#include "addToRunTimeSelectionTable.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace VF
{
defineTypeNameAndDebug(viewFactorHottel, 0);
// addToRunTimeSelectionTable(viewFactorModel, viewFactorHottel, mesh);
}
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
Foam::scalar Foam::VF::viewFactorHottel::calculateFij
(
const point& p0,
const point& p1,
const point& p2,
const point& p3
)
{
return 0.5*(mag(p2-p1) + mag(p3-p0) - mag(p2-p0) - mag(p3-p1));
}
Foam::scalarListList Foam::VF::viewFactorHottel::calculate
(
const labelListList& visibleFaceFaces,
const pointField& compactCf,
const vectorField& compactSf,
const List<List<vector>>& compactFineSf,
const List<List<point>>& compactFineCf,
const DynamicList<List<point>>& compactPoints,
const DynamicList<label>& compactPatchId
) const
{
// Fill local view factor matrix
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalarListList Fij(visibleFaceFaces.size());
forAll(visibleFaceFaces, facei)
{
if (debug > 1)
{
Pout<< "facei:" << facei << "/" << visibleFaceFaces.size()
<< endl;
}
const labelList& visibleFaces = visibleFaceFaces[facei];
Fij[facei].resize_nocopy(visibleFaces.size());
const point& dCi = compactCf[facei];
const vector& Ai = compactSf[facei];
const scalar magAi = mag(Ai);
const vector d1((Ai/magAi) ^ emptyDir_);
const vector l1(0.5*magAi/w_*d1);
const point p0(dCi + l1);
const point p1(dCi - l1);
forAll(visibleFaces, visibleFacei)
{
const label sloti = visibleFaces[visibleFacei];
const point& dCj = compactCf[sloti];
const vector& Aj = compactSf[sloti];
const scalar magAj = mag(Aj);
const vector d2((Aj/magAj) ^ emptyDir_);
const vector l2(0.5*magAj/w_*d2);
const point p2(dCj - l2);
const point p3(dCj + l2);
const scalar FijH = calculateFij(p0, p1, p2, p3);
Fij[facei][visibleFacei] = FijH/(magAi/w_);
}
}
return Fij;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::VF::viewFactorHottel::viewFactorHottel
(
const fvMesh& mesh,
const dictionary& dict
)
:
viewFactorModel(mesh, dict),
emptyDir_(vector::one),
w_(0)
{
if (mesh.nSolutionD() != 2)
{
FatalErrorInFunction
<< "Hottel crossed strings method only applicable to 2D cases"
<< exit(FatalError);
}
meshTools::constrainDirection(mesh, mesh.solutionD(), emptyDir_);
emptyDir_ = vector::one - emptyDir_;
emptyDir_.normalise();
// 2D width - assume slab
// TODO: warn wedge/axisymmetric?
w_ = mesh.bounds().span() & emptyDir_;
Info<< "\nEmpty direction: " << emptyDir_
<< "\nWidth: " << w_ << endl;
}
// ************************************************************************* //

View File

@ -0,0 +1,135 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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/>.
Class
Foam::VF::viewFactorHottel
Description
Computes view factors according to Hottel's crossed strings method.
Reference:
\verbatim
Hottel, H. C., & Saforim, A. F. (1967).
Radiative transfer.
McGraw-Hill Book Company, New York.
\endverbatim
Usage
Minimal example in \c <constant>/viewFactorsDict:
\verbatim
// Inherited entries
...
\endverbatim
The inherited entries are elaborated in:
- \link viewFactorModel.H \endlink
Note
Only applicable to 2D cases
SourceFiles
viewFactorModel.C
\*---------------------------------------------------------------------------*/
#ifndef Foam_vf_viewFactorHottel_H
#define Foam_vf_viewFactorHottel_H
#include "viewFactorModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace VF
{
/*---------------------------------------------------------------------------*\
Class viewFactorHottel Declaration
\*---------------------------------------------------------------------------*/
class viewFactorHottel
:
public viewFactorModel
{
// Private Data
//- Mesh empty direction
vector emptyDir_;
//- Mesh 2D width
scalar w_;
protected:
// Protected Member Functions
//- Calculate view factor using the double-area integral
static scalar calculateFij
(
const point& p0,
const point& p1,
const point& p2,
const point& p3
);
//- Calculate
virtual scalarListList calculate
(
const labelListList& visibleFaceFaces,
const pointField& compactCf,
const vectorField& compactSf,
const List<List<vector>>& compactFineSf,
const List<List<point>>& compactFineCf,
const DynamicList<List<point>>& compactPoints,
const DynamicList<label>& compactPatchId
) const;
public:
//- Runtime type information
TypeName("viewFactorHottel");
//- Constructor
viewFactorHottel(const fvMesh& mesh, const dictionary& dict);
//- Destructor
virtual ~viewFactorHottel() = default;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace VF
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,290 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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 "viewFactorModel.H"
#include "raySearchEngine.H"
#include "OBJstream.H"
#include "volFields.H"
#include "IOmapDistribute.H"
#include "scalarListIOList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace VF
{
defineTypeNameAndDebug(viewFactorModel, 0);
defineRunTimeSelectionTable(viewFactorModel, mesh);
}
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::VF::viewFactorModel::writeRays
(
const fileName& fName,
const pointField& compactCf,
const labelListList& visibleFaceFaces
)
{
OBJstream str(fName);
Pout<< "Writing rays to " << str.name() << endl;
forAll(visibleFaceFaces, facei)
{
const labelList& visibleSlots = visibleFaceFaces[facei];
for (const label sloti : visibleSlots)
{
str.write(linePointRef(compactCf[facei], compactCf[sloti]));
}
}
str.flush();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::VF::viewFactorModel::viewFactorModel
(
const fvMesh& mesh,
const dictionary& dict
)
:
mesh_(mesh),
searchEnginePtr_(raySearchEngine::New(mesh, dict)),
writeViewFactors_(dict.get<bool>("writeViewFactors")),
writeRays_(dict.getOrDefault<bool>("writeRays", false))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::VF::viewFactorModel::~viewFactorModel()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::VF::viewFactorModel::calculate()
{
const auto& searchEngine = *searchEnginePtr_;
labelListList visibleFaceFaces;
searchEngine.correct(visibleFaceFaces);
const auto& map = searchEngine.map();
// Construct data in compact addressing
// Compact addressing has all local data, followed by remote contributions
pointField compactCf;
vectorField compactSf;
List<List<vector>> compactFineSf;
List<List<point>> compactFineCf;
DynamicList<List<point>> compactPoints;
DynamicList<label> compactPatchId;
searchEngine.compactAddressing
(
map,
compactCf,
compactSf,
compactFineSf,
compactFineCf,
compactPoints,
compactPatchId
);
if (writeRays_)
{
// Write all rays between visible faces using .OBJ format
writeRays
(
mesh_.time().path()/"allVisibleFaces.obj",
compactCf,
visibleFaceFaces
);
}
(void)mesh_.time().cpuTimeIncrement();
Info<< "\nCalculating view factors" << endl;
scalarListIOList Fij
(
IOobject
(
"F",
mesh_.facesInstance(),
mesh_,
IOobject::NO_READ
),
calculate
(
visibleFaceFaces,
compactCf,
compactSf,
compactFineSf,
compactFineCf,
compactPoints,
compactPatchId
)
);
const label totalPatches = mesh_.boundaryMesh().nNonProcessor();
// Matrix sum in j(Fij) for each i
// Note: if enclosure sum should = 1
scalarSquareMatrix viewFactorPatch(totalPatches, Zero);
forAll(visibleFaceFaces, startFacei)
{
const scalar magAi = mag(compactSf[startFacei]);
const labelList& visibleSlots = visibleFaceFaces[startFacei];
const label patchi = compactPatchId[startFacei];
forAll(visibleSlots, visSloti)
{
const label sloti = visibleSlots[visSloti];
const label patchj = compactPatchId[sloti];
viewFactorPatch[patchi][patchj] += Fij[startFacei][visSloti]*magAi;
}
}
reduce(viewFactorPatch, sumOp<scalarSquareMatrix>());
const scalarList patchArea = searchEngine.patchAreas();
if (UPstream::master())
{
Info<< "\nPatch view factor contributions:" << nl << endl;
scalar vfSum = 0;
const auto& patchIDs = searchEngine.patchIDs();
const auto& patches = mesh_.boundaryMesh();
for (const label patchi : patchIDs)
{
Info<< " Patch " << patchi << ": " << patches[patchi].name()
<< endl;
scalar vfPatch = 0;
for (const label patchj: patchIDs)
{
scalar vf = viewFactorPatch[patchi][patchj]/patchArea[patchi];
vfPatch += vf;
Info<< " F" << patchi << patchj << ": " << vf << endl;
}
Info<< " Sum: " << vfPatch << nl << endl;
vfSum += vfPatch;
}
Info<< "Sum(all patches) = " << vfSum << endl;
}
// Write view factors matrix in list-list form
Info<< "\nWriting view factor matrix" << endl;
Fij.write();
if (writeViewFactors_)
{
Info<< "\nWriting view factor field" << endl;
auto tviewFactorField =
volScalarField::New
(
"viewFactorField",
mesh_,
dimensionedScalar(dimless, Zero)
);
auto& viewFactorField = tviewFactorField.ref();
searchEngine.interpolate(viewFactorField, Fij);
viewFactorField.write();
}
// Create globalFaceFaces needed to insert view factors
// in F to the global matrix Fmatrix
{
labelListIOList IOglobalFaceFaces
(
IOobject
(
"globalFaceFaces",
mesh_.facesInstance(),
mesh_,
IOobject::NO_READ
),
visibleFaceFaces.size()
);
forAll(IOglobalFaceFaces, facei)
{
IOglobalFaceFaces[facei] = renumber
(
searchEngine.compactToGlobal(),
visibleFaceFaces[facei]
);
}
IOglobalFaceFaces.write();
}
// Write parallel map
{
IOmapDistribute IOmapDist
(
IOobject
(
"mapDist",
mesh_.facesInstance(),
mesh_,
IOobject::NO_READ
),
std::move(map)
);
IOmapDist.write();
}
}
// ************************************************************************* //

View File

@ -0,0 +1,185 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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/>.
Namespace
Foam::VF
Description
A namespace for various \c viewFactor model implementations.
Class
Foam::VF::viewFactorModel
Description
A base class for \c viewFactor models.
Usage
Minimal example in \c <constant>/viewFactorsDict:
\verbatim
// Mandatory entries
writeViewFactors <bool>;
// Optional entries
writeRays <bool>;
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
writeViewFactors | Flag to write the view factor field | bool | yes | -
writeRays | Flag to write the ray geometry | bool | no | false
\endtable
SourceFiles
viewFactorModel.C
viewFactorModelNew.C
\*---------------------------------------------------------------------------*/
#ifndef Foam_vf_viewFactorModel_H
#define Foam_vf_viewFactorModel_H
#include "autoPtr.H"
#include "pointField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward Declarations
class fvMesh;
namespace VF
{
// Forward Declarations
class raySearchEngine;
/*---------------------------------------------------------------------------*\
Class viewFactorModel Declaration
\*---------------------------------------------------------------------------*/
class viewFactorModel
{
protected:
// Protected Data
//- Reference to the mesh database
const fvMesh& mesh_;
//- Run-time selectable ray search engine
autoPtr<raySearchEngine> searchEnginePtr_;
//- Flag to write the view factor field
bool writeViewFactors_;
//- Flag to write the ray geometry
bool writeRays_;
// Protected Member Functions
//- Write ray geometry to file
static void writeRays
(
const fileName& fName,
const pointField& compactCf,
const labelListList& visibleFaceFaces
);
//- Calculate the view factors using run-time selectable model
virtual scalarListList calculate
(
const labelListList& visibleFaceFaces,
const pointField& compactCoarseCf,
const vectorField& compactCoarseSf,
const List<List<vector>>& compactFineSf,
const List<List<point>>& compactFineCf,
const DynamicList<List<point>>& compactPoints,
const DynamicList<label>& compactPatchId
) const = 0;
public:
//- Runtime type information
TypeName("viewFactorModel");
//- Selection table
declareRunTimeSelectionTable
(
autoPtr,
viewFactorModel,
mesh,
(
const fvMesh& mesh,
const dictionary& dict
),
(mesh, dict)
);
//- Selector
static autoPtr<viewFactorModel> New
(
const fvMesh& mesh,
const dictionary& dict
);
// Generated Methods
//- No copy construct
viewFactorModel(const viewFactorModel&) = delete;
//- No copy assignment
void operator=(const viewFactorModel&) = delete;
//- Constructor
viewFactorModel(const fvMesh& mesh, const dictionary& dict);
//- Destructor
virtual ~viewFactorModel();
// Public Member Functions
//- Calculate the view factors
virtual void calculate();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace VF
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,72 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023-2024 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 "viewFactorModel.H"
#include "viewFactorHottel.H"
#include "fvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::VF::viewFactorModel> Foam::VF::viewFactorModel::New
(
const fvMesh& mesh,
const dictionary& dict
)
{
// Intercept 2D cases
if (mesh.nSolutionD() == 2)
{
Info<< "Selecting " << typeName << ": " << viewFactorHottel::typeName
<< " for 2D cases" << nl << endl;
return autoPtr<viewFactorModel>(new viewFactorHottel(mesh, dict));
}
// 3D cases - use run-time selection
const word modelType(dict.get<word>("viewFactorModel"));
Info<< "Selecting " << typeName << ": " << modelType << endl;
auto* ctorPtr = meshConstructorTable(modelType);
if (!ctorPtr)
{
FatalIOErrorInLookup
(
dict,
typeName,
modelType,
*meshConstructorTablePtr_
) << exit(FatalIOError);
}
return autoPtr<viewFactorModel>(ctorPtr(mesh, dict));
}
// ************************************************************************* //