mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-12-28 03:37:59 +00:00
- when constructing dimensioned fields that are to be zero-initialized,
it is preferrable to use a form such as
dimensionedScalar(dims, Zero)
dimensionedVector(dims, Zero)
rather than
dimensionedScalar("0", dims, 0)
dimensionedVector("zero", dims, vector::zero)
This reduces clutter and also avoids any suggestion that the name of
the dimensioned quantity has any influence on the field's name.
An even shorter version is possible. Eg,
dimensionedScalar(dims)
but reduces the clarity of meaning.
- NB: UniformDimensionedField is an exception to these style changes
since it does use the name of the dimensioned type (instead of the
regIOobject).
909 lines
25 KiB
C
909 lines
25 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
|
|
\\/ M anipulation | Copyright (C) 2016 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
|
|
viewFactorsGen
|
|
|
|
Group
|
|
grpPreProcessingUtilities
|
|
|
|
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 agglomeration.
|
|
|
|
The patches involved in the view factor calculation are taken from the
|
|
boundary file and should be part on the group viewFactorWall. ie.:
|
|
|
|
floor
|
|
{
|
|
type wall;
|
|
inGroups 2(wall viewFactorWall);
|
|
nFaces 100;
|
|
startFace 3100;
|
|
}
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "argList.H"
|
|
#include "fvMesh.H"
|
|
#include "volFields.H"
|
|
#include "surfaceFields.H"
|
|
#include "distributedTriSurfaceMesh.H"
|
|
#include "meshTools.H"
|
|
|
|
#include "uindirectPrimitivePatch.H"
|
|
#include "DynamicField.H"
|
|
#include "unitConversion.H"
|
|
|
|
#include "scalarMatrices.H"
|
|
#include "labelListIOList.H"
|
|
#include "scalarListIOList.H"
|
|
|
|
#include "singleCellFvMesh.H"
|
|
|
|
#include "IOmapDistribute.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++;
|
|
}
|
|
|
|
//striSurfaceToAgglom.resize(localTriFaceI-1);
|
|
|
|
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;
|
|
}
|
|
}
|
|
str.flush();
|
|
|
|
Pout<< "cmd: objToVTK " << fName.c_str() << endl;
|
|
|
|
stringList cmd{"objToVTK", fName, fName.lessExt().ext("vtk")};
|
|
Foam::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 word viewFactorWall("viewFactorWall");
|
|
|
|
const bool writeViewFactors =
|
|
viewFactorDict.lookupOrDefault("writeViewFactorMatrix", false);
|
|
|
|
const bool dumpRays =
|
|
viewFactorDict.lookupOrDefault("dumpRays", false);
|
|
|
|
const label debug = viewFactorDict.lookupOrDefault<label>("debug", 0);
|
|
|
|
// 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.findIndices(viewFactorWall));
|
|
forAll(viewFactorsPatches, i)
|
|
{
|
|
label patchI = viewFactorsPatches[i];
|
|
nCoarseFaces += coarsePatches[patchI].size();
|
|
nFineFaces += patches[patchI].size();
|
|
}
|
|
|
|
// 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);
|
|
labelHashSet includePatches;
|
|
|
|
forAll(viewFactorsPatches, i)
|
|
{
|
|
const label patchID = viewFactorsPatches[i];
|
|
|
|
const polyPatch& pp = patches[patchID];
|
|
const labelList& agglom = finalAgglom[patchID];
|
|
|
|
includePatches.insert(patchID);
|
|
|
|
if (agglom.size() > 0)
|
|
{
|
|
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];
|
|
|
|
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 and rayEndFace 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);
|
|
|
|
// 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];
|
|
if (agglom.size() > 0)
|
|
{
|
|
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(dimless, Zero)
|
|
);
|
|
|
|
label compactI = 0;
|
|
|
|
volScalarField::Boundary& vfbf = viewFactorField.boundaryFieldRef();
|
|
forAll(viewFactorsPatches, i)
|
|
{
|
|
label patchID = viewFactorsPatches[i];
|
|
const labelList& agglom = finalAgglom[patchID];
|
|
if (agglom.size() > 0)
|
|
{
|
|
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];
|
|
vfbf[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()
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
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
|
|
),
|
|
std::move(globalFaceFaces)
|
|
);
|
|
IOglobalFaceFaces.write();
|
|
|
|
|
|
labelListIOList IOvisibleFaceFaces
|
|
(
|
|
IOobject
|
|
(
|
|
"visibleFaceFaces",
|
|
mesh.facesInstance(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
std::move(visibleFaceFaces)
|
|
);
|
|
IOvisibleFaceFaces.write();
|
|
|
|
IOmapDistribute IOmapDist
|
|
(
|
|
IOobject
|
|
(
|
|
"mapDist",
|
|
mesh.facesInstance(),
|
|
mesh,
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
std::move(map)
|
|
);
|
|
|
|
IOmapDist.write();
|
|
|
|
Info<< "End\n" << endl;
|
|
return 0;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|