955 lines
26 KiB
C
955 lines
26 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM is free software: you can redistribute it and/or modify it
|
|
under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
Application
|
|
viewFactorsGen
|
|
|
|
Description
|
|
View factors are calculated based on a face agglomeration array
|
|
(finalAgglom generated by faceAgglomerate utility).
|
|
|
|
Each view factor between the agglomerated faces i and j (Fij) is calculated
|
|
using a double integral of the sub-areas composing the agglomaration.
|
|
|
|
The patches involved in the view factor calculation are taken from the Qr
|
|
volScalarField (radiative flux) when is greyDiffusiveRadiationViewFactor
|
|
otherwise they are not included.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "argList.H"
|
|
#include "Time.H"
|
|
#include "fvMesh.H"
|
|
#include "singleCellFvMesh.H"
|
|
#include "volFields.H"
|
|
#include "surfaceFields.H"
|
|
#include "fixedValueFvPatchFields.H"
|
|
#include "distributedTriSurfaceMesh.H"
|
|
#include "cyclicAMIPolyPatch.H"
|
|
#include "mapDistribute.H"
|
|
#include "meshTools.H"
|
|
#include "uindirectPrimitivePatch.H"
|
|
#include "DynamicField.H"
|
|
#include "scalarMatrices.H"
|
|
#include "scalarListIOList.H"
|
|
|
|
using namespace Foam;
|
|
|
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
|
|
|
triSurface triangulate
|
|
(
|
|
const polyBoundaryMesh& bMesh,
|
|
const labelHashSet& includePatches,
|
|
const labelListIOList& finalAgglom,
|
|
labelList& triSurfaceToAgglom,
|
|
const globalIndex& globalNumbering,
|
|
const polyBoundaryMesh& coarsePatches
|
|
)
|
|
{
|
|
const polyMesh& mesh = bMesh.mesh();
|
|
|
|
// Storage for surfaceMesh. Size estimate.
|
|
DynamicList<labelledTri> triangles
|
|
(
|
|
mesh.nFaces() - mesh.nInternalFaces()
|
|
);
|
|
|
|
label newPatchI = 0;
|
|
label localTriFaceI = 0;
|
|
|
|
forAllConstIter(labelHashSet, includePatches, iter)
|
|
{
|
|
const label patchI = iter.key();
|
|
const polyPatch& patch = bMesh[patchI];
|
|
const pointField& points = patch.points();
|
|
|
|
label nTriTotal = 0;
|
|
|
|
forAll(patch, patchFaceI)
|
|
{
|
|
const face& f = patch[patchFaceI];
|
|
|
|
faceList triFaces(f.nTriangles(points));
|
|
|
|
label nTri = 0;
|
|
|
|
f.triangles(points, nTri, triFaces);
|
|
|
|
forAll(triFaces, triFaceI)
|
|
{
|
|
const face& f = triFaces[triFaceI];
|
|
|
|
triangles.append(labelledTri(f[0], f[1], f[2], newPatchI));
|
|
|
|
nTriTotal++;
|
|
|
|
triSurfaceToAgglom[localTriFaceI++] = globalNumbering.toGlobal
|
|
(
|
|
Pstream::myProcNo(),
|
|
finalAgglom[patchI][patchFaceI]
|
|
+ coarsePatches[patchI].start()
|
|
);
|
|
}
|
|
}
|
|
|
|
newPatchI++;
|
|
}
|
|
|
|
triSurfaceToAgglom.resize(localTriFaceI);
|
|
|
|
triangles.shrink();
|
|
|
|
// Create globally numbered tri surface
|
|
triSurface rawSurface(triangles, mesh.points());
|
|
|
|
// Create locally numbered tri surface
|
|
triSurface surface
|
|
(
|
|
rawSurface.localFaces(),
|
|
rawSurface.localPoints()
|
|
);
|
|
|
|
// Add patch names to surface
|
|
surface.patches().setSize(newPatchI);
|
|
|
|
newPatchI = 0;
|
|
|
|
forAllConstIter(labelHashSet, includePatches, iter)
|
|
{
|
|
const label patchI = iter.key();
|
|
const polyPatch& patch = bMesh[patchI];
|
|
|
|
surface.patches()[newPatchI].index() = patchI;
|
|
surface.patches()[newPatchI].name() = patch.name();
|
|
surface.patches()[newPatchI].geometricType() = patch.type();
|
|
|
|
newPatchI++;
|
|
}
|
|
|
|
return surface;
|
|
}
|
|
|
|
|
|
void writeRays
|
|
(
|
|
const fileName& fName,
|
|
const pointField& compactCf,
|
|
const pointField& myFc,
|
|
const labelListList& visibleFaceFaces
|
|
)
|
|
{
|
|
OFstream str(fName);
|
|
label vertI = 0;
|
|
|
|
Pout<< "Dumping rays to " << str.name() << endl;
|
|
|
|
forAll(myFc, faceI)
|
|
{
|
|
const labelList visFaces = visibleFaceFaces[faceI];
|
|
forAll(visFaces, faceRemote)
|
|
{
|
|
label compactI = visFaces[faceRemote];
|
|
const point& remoteFc = compactCf[compactI];
|
|
|
|
meshTools::writeOBJ(str, myFc[faceI]);
|
|
vertI++;
|
|
meshTools::writeOBJ(str, remoteFc);
|
|
vertI++;
|
|
str << "l " << vertI-1 << ' ' << vertI << nl;
|
|
}
|
|
}
|
|
string cmd("objToVTK " + fName + " " + fName.lessExt() + ".vtk");
|
|
Pout<< "cmd:" << cmd << endl;
|
|
system(cmd);
|
|
}
|
|
|
|
|
|
scalar calculateViewFactorFij
|
|
(
|
|
const vector& i,
|
|
const vector& j,
|
|
const vector& dAi,
|
|
const vector& dAj
|
|
)
|
|
{
|
|
vector r = i - j;
|
|
scalar rMag = mag(r);
|
|
|
|
if (rMag > SMALL)
|
|
{
|
|
scalar dAiMag = mag(dAi);
|
|
scalar dAjMag = mag(dAj);
|
|
|
|
vector ni = dAi/dAiMag;
|
|
vector nj = dAj/dAjMag;
|
|
scalar cosThetaJ = mag(nj & r)/rMag;
|
|
scalar cosThetaI = mag(ni & r)/rMag;
|
|
|
|
return
|
|
(
|
|
(cosThetaI*cosThetaJ*dAjMag*dAiMag)
|
|
/(sqr(rMag)*constant::mathematical::pi)
|
|
);
|
|
}
|
|
else
|
|
{
|
|
return 0;
|
|
}
|
|
}
|
|
|
|
|
|
void insertMatrixElements
|
|
(
|
|
const globalIndex& globalNumbering,
|
|
const label fromProcI,
|
|
const labelListList& globalFaceFaces,
|
|
const scalarListList& viewFactors,
|
|
scalarSquareMatrix& matrix
|
|
)
|
|
{
|
|
forAll(viewFactors, faceI)
|
|
{
|
|
const scalarList& vf = viewFactors[faceI];
|
|
const labelList& globalFaces = globalFaceFaces[faceI];
|
|
|
|
label globalI = globalNumbering.toGlobal(fromProcI, faceI);
|
|
forAll(globalFaces, i)
|
|
{
|
|
matrix[globalI][globalFaces[i]] = vf[i];
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
|
|
|
int main(int argc, char *argv[])
|
|
{
|
|
#include "addRegionOption.H"
|
|
#include "setRootCase.H"
|
|
#include "createTime.H"
|
|
#include "createNamedMesh.H"
|
|
|
|
// Read view factor dictionary
|
|
IOdictionary viewFactorDict
|
|
(
|
|
IOobject
|
|
(
|
|
"viewFactorsDict",
|
|
runTime.constant(),
|
|
mesh,
|
|
IOobject::MUST_READ_IF_MODIFIED,
|
|
IOobject::NO_WRITE
|
|
)
|
|
);
|
|
|
|
const bool writeViewFactors =
|
|
viewFactorDict.lookupOrDefault<bool>("writeViewFactorMatrix", false);
|
|
|
|
const bool dumpRays =
|
|
viewFactorDict.lookupOrDefault<bool>("dumpRays", false);
|
|
|
|
const label debug = viewFactorDict.lookupOrDefault<label>("debug", 0);
|
|
|
|
volScalarField Qr
|
|
(
|
|
IOobject
|
|
(
|
|
"Qr",
|
|
runTime.timeName(),
|
|
mesh,
|
|
IOobject::MUST_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh
|
|
);
|
|
|
|
// Read agglomeration map
|
|
labelListIOList finalAgglom
|
|
(
|
|
IOobject
|
|
(
|
|
"finalAgglom",
|
|
mesh.facesInstance(),
|
|
mesh,
|
|
IOobject::MUST_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
)
|
|
);
|
|
|
|
// Create the coarse mesh using agglomeration
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
if (debug)
|
|
{
|
|
Pout << "\nCreating single cell mesh..." << endl;
|
|
}
|
|
|
|
singleCellFvMesh coarseMesh
|
|
(
|
|
IOobject
|
|
(
|
|
"coarse:" + mesh.name(),
|
|
runTime.timeName(),
|
|
runTime,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh,
|
|
finalAgglom
|
|
);
|
|
|
|
if (debug)
|
|
{
|
|
Pout << "\nCreated single cell mesh..." << endl;
|
|
}
|
|
|
|
|
|
// Calculate total number of fine and coarse faces
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
label nCoarseFaces = 0; //total number of coarse faces
|
|
label nFineFaces = 0; //total number of fine faces
|
|
|
|
const polyBoundaryMesh& patches = mesh.boundaryMesh();
|
|
const polyBoundaryMesh& coarsePatches = coarseMesh.boundaryMesh();
|
|
|
|
labelList viewFactorsPatches(patches.size());
|
|
|
|
const volScalarField::GeometricBoundaryField& Qrb = Qr.boundaryField();
|
|
|
|
label count = 0;
|
|
forAll(Qrb, patchI)
|
|
{
|
|
const polyPatch& pp = patches[patchI];
|
|
const fvPatchScalarField& QrpI = Qrb[patchI];
|
|
|
|
if ((isA<fixedValueFvPatchScalarField>(QrpI)) && (pp.size() > 0))
|
|
{
|
|
viewFactorsPatches[count] = QrpI.patch().index();
|
|
nCoarseFaces += coarsePatches[patchI].size();
|
|
nFineFaces += patches[patchI].size();
|
|
count ++;
|
|
}
|
|
}
|
|
|
|
viewFactorsPatches.resize(count);
|
|
|
|
// total number of coarse faces
|
|
label totalNCoarseFaces = nCoarseFaces;
|
|
|
|
reduce(totalNCoarseFaces, sumOp<label>());
|
|
|
|
if (Pstream::master())
|
|
{
|
|
Info << "\nTotal number of coarse faces: "<< totalNCoarseFaces << endl;
|
|
}
|
|
|
|
if (Pstream::master() && debug)
|
|
{
|
|
Pout << "\nView factor patches included in the calculation : "
|
|
<< viewFactorsPatches << endl;
|
|
}
|
|
|
|
// Collect local Cf and Sf on coarse mesh
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
DynamicList<point> localCoarseCf(nCoarseFaces);
|
|
DynamicList<point> localCoarseSf(nCoarseFaces);
|
|
|
|
DynamicList<label> localAgg(nCoarseFaces);
|
|
|
|
forAll(viewFactorsPatches, i)
|
|
{
|
|
const label patchID = viewFactorsPatches[i];
|
|
|
|
const polyPatch& pp = patches[patchID];
|
|
const labelList& agglom = finalAgglom[patchID];
|
|
label nAgglom = max(agglom)+1;
|
|
labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
|
|
const labelList& coarsePatchFace = coarseMesh.patchFaceMap()[patchID];
|
|
|
|
const pointField& coarseCf = coarseMesh.Cf().boundaryField()[patchID];
|
|
const pointField& coarseSf = coarseMesh.Sf().boundaryField()[patchID];
|
|
|
|
labelHashSet includePatches;
|
|
includePatches.insert(patchID);
|
|
|
|
forAll(coarseCf, faceI)
|
|
{
|
|
point cf = coarseCf[faceI];
|
|
|
|
const label coarseFaceI = coarsePatchFace[faceI];
|
|
const labelList& fineFaces = coarseToFine[coarseFaceI];
|
|
const label agglomI =
|
|
agglom[fineFaces[0]] + coarsePatches[patchID].start();
|
|
|
|
// Construct single face
|
|
uindirectPrimitivePatch upp
|
|
(
|
|
UIndirectList<face>(pp, fineFaces),
|
|
pp.points()
|
|
);
|
|
|
|
|
|
List<point> availablePoints
|
|
(
|
|
upp.faceCentres().size()
|
|
+ upp.localPoints().size()
|
|
);
|
|
|
|
SubList<point>
|
|
(
|
|
availablePoints,
|
|
upp.faceCentres().size()
|
|
) = upp.faceCentres();
|
|
|
|
SubList<point>
|
|
(
|
|
availablePoints,
|
|
upp.localPoints().size(),
|
|
upp.faceCentres().size()
|
|
) = upp.localPoints();
|
|
|
|
point cfo = cf;
|
|
scalar dist = GREAT;
|
|
forAll(availablePoints, iPoint)
|
|
{
|
|
point cfFine = availablePoints[iPoint];
|
|
if (mag(cfFine-cfo) < dist)
|
|
{
|
|
dist = mag(cfFine-cfo);
|
|
cf = cfFine;
|
|
}
|
|
}
|
|
|
|
point sf = coarseSf[faceI];
|
|
localCoarseCf.append(cf);
|
|
localCoarseSf.append(sf);
|
|
localAgg.append(agglomI);
|
|
}
|
|
}
|
|
|
|
|
|
// Distribute local coarse Cf and Sf for shooting rays
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
List<pointField> remoteCoarseCf(Pstream::nProcs());
|
|
List<pointField> remoteCoarseSf(Pstream::nProcs());
|
|
List<labelField> remoteCoarseAgg(Pstream::nProcs());
|
|
|
|
remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
|
|
remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
|
|
remoteCoarseAgg[Pstream::myProcNo()] = localAgg;
|
|
|
|
Pstream::gatherList(remoteCoarseCf);
|
|
Pstream::scatterList(remoteCoarseCf);
|
|
Pstream::gatherList(remoteCoarseSf);
|
|
Pstream::scatterList(remoteCoarseSf);
|
|
Pstream::gatherList(remoteCoarseAgg);
|
|
Pstream::scatterList(remoteCoarseAgg);
|
|
|
|
|
|
globalIndex globalNumbering(nCoarseFaces);
|
|
|
|
// Set up searching engine for obstacles
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
#include "searchingEngine.H"
|
|
|
|
|
|
// Determine rays between coarse face centres
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
DynamicList<label> rayStartFace(nCoarseFaces + 0.01*nCoarseFaces);
|
|
|
|
DynamicList<label> rayEndFace(rayStartFace.size());
|
|
|
|
|
|
// Return rayStartFace in local index andrayEndFace in global index
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
#include "shootRays.H"
|
|
|
|
// Calculate number of visible faces from local index
|
|
labelList nVisibleFaceFaces(nCoarseFaces, 0);
|
|
|
|
forAll(rayStartFace, i)
|
|
{
|
|
nVisibleFaceFaces[rayStartFace[i]]++;
|
|
}
|
|
|
|
labelListList visibleFaceFaces(nCoarseFaces);
|
|
|
|
label nViewFactors = 0;
|
|
forAll(nVisibleFaceFaces, faceI)
|
|
{
|
|
visibleFaceFaces[faceI].setSize(nVisibleFaceFaces[faceI]);
|
|
nViewFactors += nVisibleFaceFaces[faceI];
|
|
}
|
|
|
|
|
|
// - 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
|
|
|
|
List<Map<label>> compactMap(Pstream::nProcs());
|
|
|
|
mapDistribute map(globalNumbering, rayEndFace, compactMap);
|
|
|
|
labelListIOList IOsubMap
|
|
(
|
|
IOobject
|
|
(
|
|
"subMap",
|
|
mesh.facesInstance(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
map.subMap()
|
|
);
|
|
IOsubMap.write();
|
|
|
|
|
|
labelListIOList IOconstructMap
|
|
(
|
|
IOobject
|
|
(
|
|
"constructMap",
|
|
mesh.facesInstance(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
map.constructMap()
|
|
);
|
|
IOconstructMap.write();
|
|
|
|
|
|
IOList<label> consMapDim
|
|
(
|
|
IOobject
|
|
(
|
|
"constructMapDim",
|
|
mesh.facesInstance(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
List<label>(1, map.constructSize())
|
|
);
|
|
consMapDim.write();
|
|
|
|
|
|
// visibleFaceFaces has:
|
|
// (local face, local viewed face) = compact viewed face
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
nVisibleFaceFaces = 0;
|
|
forAll(rayStartFace, i)
|
|
{
|
|
label faceI = rayStartFace[i];
|
|
label compactI = rayEndFace[i];
|
|
visibleFaceFaces[faceI][nVisibleFaceFaces[faceI]++] = compactI;
|
|
}
|
|
|
|
|
|
// Construct data in compact addressing
|
|
// I need coarse Sf (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
pointField compactCoarseCf(map.constructSize(), Zero);
|
|
pointField compactCoarseSf(map.constructSize(), Zero);
|
|
List<List<point>> compactFineSf(map.constructSize());
|
|
List<List<point>> compactFineCf(map.constructSize());
|
|
|
|
DynamicList<label> compactPatchId(map.constructSize());
|
|
|
|
// Insert my coarse local values
|
|
SubList<point>(compactCoarseSf, nCoarseFaces) = localCoarseSf;
|
|
SubList<point>(compactCoarseCf, nCoarseFaces) = localCoarseCf;
|
|
|
|
// Insert my fine local values
|
|
label compactI = 0;
|
|
forAll(viewFactorsPatches, i)
|
|
{
|
|
label patchID = viewFactorsPatches[i];
|
|
const labelList& agglom = finalAgglom[patchID];
|
|
label nAgglom = max(agglom)+1;
|
|
labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
|
|
const labelList& coarsePatchFace = coarseMesh.patchFaceMap()[patchID];
|
|
|
|
forAll(coarseToFine, coarseI)
|
|
{
|
|
compactPatchId.append(patchID);
|
|
List<point>& fineCf = compactFineCf[compactI];
|
|
List<point>& fineSf = compactFineSf[compactI++];
|
|
|
|
const label coarseFaceI = coarsePatchFace[coarseI];
|
|
const labelList& fineFaces = coarseToFine[coarseFaceI];
|
|
|
|
fineCf.setSize(fineFaces.size());
|
|
fineSf.setSize(fineFaces.size());
|
|
|
|
fineCf = UIndirectList<point>
|
|
(
|
|
mesh.Cf().boundaryField()[patchID],
|
|
coarseToFine[coarseFaceI]
|
|
);
|
|
fineSf = UIndirectList<point>
|
|
(
|
|
mesh.Sf().boundaryField()[patchID],
|
|
coarseToFine[coarseFaceI]
|
|
);
|
|
}
|
|
}
|
|
|
|
// Do all swapping
|
|
map.distribute(compactCoarseSf);
|
|
map.distribute(compactCoarseCf);
|
|
map.distribute(compactFineCf);
|
|
map.distribute(compactFineSf);
|
|
|
|
map.distribute(compactPatchId);
|
|
|
|
|
|
// Plot all rays between visible faces.
|
|
if (dumpRays)
|
|
{
|
|
writeRays
|
|
(
|
|
runTime.path()/"allVisibleFaces.obj",
|
|
compactCoarseCf,
|
|
remoteCoarseCf[Pstream::myProcNo()],
|
|
visibleFaceFaces
|
|
);
|
|
}
|
|
|
|
|
|
// Fill local view factor matrix
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
scalarListIOList F
|
|
(
|
|
IOobject
|
|
(
|
|
"F",
|
|
mesh.facesInstance(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
nCoarseFaces
|
|
);
|
|
|
|
label totalPatches = coarsePatches.size();
|
|
reduce(totalPatches, maxOp<label>());
|
|
|
|
// Matrix sum in j(Fij) for each i (if enclosure sum = 1)
|
|
scalarSquareMatrix sumViewFactorPatch
|
|
(
|
|
totalPatches,
|
|
0.0
|
|
);
|
|
|
|
scalarList patchArea(totalPatches, 0.0);
|
|
|
|
if (Pstream::master())
|
|
{
|
|
Info<< "\nCalculating view factors..." << endl;
|
|
}
|
|
|
|
if (mesh.nSolutionD() == 3)
|
|
{
|
|
forAll(localCoarseSf, coarseFaceI)
|
|
{
|
|
const List<point>& localFineSf = compactFineSf[coarseFaceI];
|
|
const vector Ai = sum(localFineSf);
|
|
const List<point>& localFineCf = compactFineCf[coarseFaceI];
|
|
const label fromPatchId = compactPatchId[coarseFaceI];
|
|
patchArea[fromPatchId] += mag(Ai);
|
|
|
|
const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
|
|
|
|
forAll(visCoarseFaces, visCoarseFaceI)
|
|
{
|
|
F[coarseFaceI].setSize(visCoarseFaces.size());
|
|
label compactJ = visCoarseFaces[visCoarseFaceI];
|
|
const List<point>& remoteFineSj = compactFineSf[compactJ];
|
|
const List<point>& remoteFineCj = compactFineCf[compactJ];
|
|
|
|
const label toPatchId = compactPatchId[compactJ];
|
|
|
|
scalar Fij = 0;
|
|
forAll(localFineSf, i)
|
|
{
|
|
const vector& dAi = localFineSf[i];
|
|
const vector& dCi = localFineCf[i];
|
|
|
|
forAll(remoteFineSj, j)
|
|
{
|
|
const vector& dAj = remoteFineSj[j];
|
|
const vector& dCj = remoteFineCj[j];
|
|
|
|
scalar dIntFij = calculateViewFactorFij
|
|
(
|
|
dCi,
|
|
dCj,
|
|
dAi,
|
|
dAj
|
|
);
|
|
|
|
Fij += dIntFij;
|
|
}
|
|
}
|
|
F[coarseFaceI][visCoarseFaceI] = Fij/mag(Ai);
|
|
sumViewFactorPatch[fromPatchId][toPatchId] += Fij;
|
|
}
|
|
}
|
|
}
|
|
else if (mesh.nSolutionD() == 2)
|
|
{
|
|
const boundBox& box = mesh.bounds();
|
|
const Vector<label>& dirs = mesh.geometricD();
|
|
vector emptyDir = Zero;
|
|
forAll(dirs, i)
|
|
{
|
|
if (dirs[i] == -1)
|
|
{
|
|
emptyDir[i] = 1.0;
|
|
}
|
|
}
|
|
|
|
scalar wideBy2 = (box.span() & emptyDir)*2.0;
|
|
|
|
forAll(localCoarseSf, coarseFaceI)
|
|
{
|
|
const vector& Ai = localCoarseSf[coarseFaceI];
|
|
const vector& Ci = localCoarseCf[coarseFaceI];
|
|
vector Ain = Ai/mag(Ai);
|
|
vector R1i = Ci + (mag(Ai)/wideBy2)*(Ain ^ emptyDir);
|
|
vector R2i = Ci - (mag(Ai)/wideBy2)*(Ain ^ emptyDir) ;
|
|
|
|
const label fromPatchId = compactPatchId[coarseFaceI];
|
|
patchArea[fromPatchId] += mag(Ai);
|
|
|
|
const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
|
|
forAll(visCoarseFaces, visCoarseFaceI)
|
|
{
|
|
F[coarseFaceI].setSize(visCoarseFaces.size());
|
|
label compactJ = visCoarseFaces[visCoarseFaceI];
|
|
const vector& Aj = compactCoarseSf[compactJ];
|
|
const vector& Cj = compactCoarseCf[compactJ];
|
|
|
|
const label toPatchId = compactPatchId[compactJ];
|
|
|
|
vector Ajn = Aj/mag(Aj);
|
|
vector R1j = Cj + (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
|
|
vector R2j = Cj - (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
|
|
|
|
scalar d1 = mag(R1i - R2j);
|
|
scalar d2 = mag(R2i - R1j);
|
|
scalar s1 = mag(R1i - R1j);
|
|
scalar s2 = mag(R2i - R2j);
|
|
|
|
scalar Fij = mag((d1 + d2) - (s1 + s2))/(4.0*mag(Ai)/wideBy2);
|
|
|
|
F[coarseFaceI][visCoarseFaceI] = Fij;
|
|
sumViewFactorPatch[fromPatchId][toPatchId] += Fij*mag(Ai);
|
|
}
|
|
}
|
|
}
|
|
|
|
if (Pstream::master())
|
|
{
|
|
Info << "Writing view factor matrix..." << endl;
|
|
}
|
|
|
|
// Write view factors matrix in listlist form
|
|
F.write();
|
|
|
|
reduce(sumViewFactorPatch, sumOp<scalarSquareMatrix>());
|
|
reduce(patchArea, sumOp<scalarList>());
|
|
|
|
|
|
if (Pstream::master() && debug)
|
|
{
|
|
forAll(viewFactorsPatches, i)
|
|
{
|
|
label patchI = viewFactorsPatches[i];
|
|
forAll(viewFactorsPatches, i)
|
|
{
|
|
label patchJ = viewFactorsPatches[i];
|
|
Info << "F" << patchI << patchJ << ": "
|
|
<< sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
|
|
<< endl;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
if (writeViewFactors)
|
|
{
|
|
volScalarField viewFactorField
|
|
(
|
|
IOobject
|
|
(
|
|
"viewFactorField",
|
|
mesh.time().timeName(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh,
|
|
dimensionedScalar("viewFactorField", dimless, 0)
|
|
);
|
|
|
|
label compactI = 0;
|
|
forAll(viewFactorsPatches, i)
|
|
{
|
|
label patchID = viewFactorsPatches[i];
|
|
const labelList& agglom = finalAgglom[patchID];
|
|
label nAgglom = max(agglom)+1;
|
|
labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
|
|
const labelList& coarsePatchFace =
|
|
coarseMesh.patchFaceMap()[patchID];
|
|
|
|
forAll(coarseToFine, coarseI)
|
|
{
|
|
const scalar Fij = sum(F[compactI]);
|
|
const label coarseFaceID = coarsePatchFace[coarseI];
|
|
const labelList& fineFaces = coarseToFine[coarseFaceID];
|
|
forAll(fineFaces, fineId)
|
|
{
|
|
const label faceID = fineFaces[fineId];
|
|
viewFactorField.boundaryField()[patchID][faceID] = Fij;
|
|
}
|
|
compactI++;
|
|
}
|
|
}
|
|
viewFactorField.write();
|
|
}
|
|
|
|
|
|
// Invert compactMap (from processor+localface to compact) to go
|
|
// from compact to processor+localface (expressed as a globalIndex)
|
|
// globalIndex globalCoarFaceNum(coarseMesh.nFaces());
|
|
labelList compactToGlobal(map.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];
|
|
|
|
forAllConstIter(Map<label>, localToCompactMap, iter)
|
|
{
|
|
compactToGlobal[iter()] = globalNumbering.toGlobal
|
|
(
|
|
procI,
|
|
iter.key()
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
if (Pstream::master())
|
|
{
|
|
scalarSquareMatrix Fmatrix(totalNCoarseFaces, 0.0);
|
|
|
|
labelListList globalFaceFaces(visibleFaceFaces.size());
|
|
|
|
// Create globalFaceFaces needed to insert view factors
|
|
// in F to the global matrix Fmatrix
|
|
forAll(globalFaceFaces, faceI)
|
|
{
|
|
globalFaceFaces[faceI] = renumber
|
|
(
|
|
compactToGlobal,
|
|
visibleFaceFaces[faceI]
|
|
);
|
|
}
|
|
|
|
labelListIOList IOglobalFaceFaces
|
|
(
|
|
IOobject
|
|
(
|
|
"globalFaceFaces",
|
|
mesh.facesInstance(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
globalFaceFaces
|
|
);
|
|
IOglobalFaceFaces.write();
|
|
}
|
|
else
|
|
{
|
|
labelListList globalFaceFaces(visibleFaceFaces.size());
|
|
forAll(globalFaceFaces, faceI)
|
|
{
|
|
globalFaceFaces[faceI] = renumber
|
|
(
|
|
compactToGlobal,
|
|
visibleFaceFaces[faceI]
|
|
);
|
|
}
|
|
|
|
labelListIOList IOglobalFaceFaces
|
|
(
|
|
IOobject
|
|
(
|
|
"globalFaceFaces",
|
|
mesh.facesInstance(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
globalFaceFaces
|
|
);
|
|
|
|
IOglobalFaceFaces.write();
|
|
}
|
|
|
|
Info<< "End\n" << endl;
|
|
return 0;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|