Merge remote branch 'OpenCFD/master' into olesenm

Conflicts:
	bin/foamCleanPath
This commit is contained in:
Mark Olesen
2011-02-22 19:24:31 +01:00
156 changed files with 8498 additions and 919 deletions

View File

@ -60,6 +60,8 @@ Usage
\param -literalRE \n
Do not interpret regular expressions; treat them as any other keyword.
\param -enableFunctionEntries \n
By default all dictionary preprocessing of fields is disabled
\*---------------------------------------------------------------------------*/
@ -255,6 +257,11 @@ int main(int argc, char *argv[])
"literalRE",
"treat regular expressions literally (ie, as a keyword)"
);
argList::addBoolOption
(
"enableFunctionEntries",
"enable expansion of dictionary directives - #include, #codeStream etc"
);
#include "addRegionOption.H"
#include "setRootCase.H"
@ -271,6 +278,20 @@ int main(int argc, char *argv[])
<< " not present." << endl;
}
const bool enableEntries = args.optionFound("enableFunctionEntries");
if (enableEntries)
{
Info<< "Allowing dictionary preprocessing ('#include', '#codeStream')."
<< endl;
}
int oldFlag = entry::disableFunctionEntries;
if (!enableEntries)
{
// By default disable dictionary expansion for fields
entry::disableFunctionEntries = 1;
}
fileName regionPrefix = "";
if (regionName != fvMesh::defaultRegion)
@ -410,6 +431,7 @@ int main(int argc, char *argv[])
false
)
);
const_cast<word&>(IOdictionary::typeName) = oldTypeName;
// Fake type back to what was in field
const_cast<word&>(fieldDict.type()) = fieldDict.headerClassName();
@ -429,6 +451,8 @@ int main(int argc, char *argv[])
}
}
entry::disableFunctionEntries = oldFlag;
Info<< endl;
Info<< "End\n" << endl;

View File

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

View File

@ -0,0 +1,11 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/OpenFOAM/lnInclude \
-I$(LIB_SRC)/fvAgglomerationMethods/pairPatchAgglomeration/lnInclude
EXE_LIBS = \
-lOpenFOAM \
-lfiniteVolume \
-lpairPatchAgglomeration \
-ltriSurface \
-lmeshTools

View File

@ -0,0 +1,192 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ 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
faceAgglomerate
Description
Agglomerate boundary faces using the pairPatchAgglomeration algorithm.
It writes a map from the fine to coarse grid.
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "fvMesh.H"
#include "Time.H"
#include "volFields.H"
#include "CompactListList.H"
#include "unitConversion.H"
#include "pairPatchAgglomeration.H"
#include "labelListIOList.H"
#include "syncTools.H"
using namespace Foam;
// Main program:
int main(int argc, char *argv[])
{
#include "addRegionOption.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createNamedMesh.H"
const polyBoundaryMesh& patches = mesh.boundaryMesh();
labelListIOList finalAgglom
(
IOobject
(
"finalAgglom",
mesh.facesInstance(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
patches.size()
);
// Read view factor dictionary
IOdictionary viewFactorDict
(
IOobject
(
"viewFactorsDict",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
bool writeAgglo =
readBool(viewFactorDict.lookup("writeFacesAgglomeration"));
const polyBoundaryMesh& boundary = mesh.boundaryMesh();
forAll(boundary, patchId)
{
const polyPatch& pp = boundary[patchId];
label patchI = pp.index();
finalAgglom[patchI].setSize(pp.size(), 0);
if (pp.size() > 0 && !pp.coupled())
{
if (viewFactorDict.found(pp.name()))
{
Info << "\nAgglomerating name : " << pp.name() << endl;
pairPatchAgglomeration agglomObject
(
pp,
viewFactorDict.subDict(pp.name())
);
agglomObject.agglomerate();
finalAgglom[patchI] =
agglomObject.restrictTopBottomAddressing();
}
else
{
FatalErrorIn
(
"main(int argc, char *argv[])"
) << pp.name()
<< " not found in dictionary : "
<< viewFactorDict.name()
<< exit(FatalError);
}
}
}
// Sync agglomeration across coupled patches
labelList nbrAgglom(mesh.nFaces()-mesh.nInternalFaces(), -1);
forAll(boundary, patchId)
{
const polyPatch& pp = boundary[patchId];
if (pp.coupled())
{
finalAgglom[patchId] = identity(pp.size());
forAll(pp, i)
{
nbrAgglom[pp.start()-mesh.nInternalFaces()+i] =
finalAgglom[patchId][i];
}
}
}
syncTools::swapBoundaryFaceList(mesh, nbrAgglom);
forAll(boundary, patchId)
{
const polyPatch& pp = boundary[patchId];
if (pp.coupled() && !refCast<const coupledPolyPatch>(pp).owner())
{
forAll(pp, i)
{
finalAgglom[patchId][i] =
nbrAgglom[pp.start()-mesh.nInternalFaces()+i];
}
}
}
finalAgglom.write();
if (writeAgglo)
{
volScalarField facesAgglomeration
(
IOobject
(
"facesAgglomeration",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("facesAgglomeration", dimless, 0)
);
forAll(boundary, patchId)
{
fvPatchScalarField& bFacesAgglomeration =
facesAgglomeration.boundaryField()[patchId];
forAll(bFacesAgglomeration, j)
{
bFacesAgglomeration[j] = finalAgglom[patchId][j];
}
}
Info << "\nWriting facesAgglomeration..." << endl;
facesAgglomeration.write();
}
Info<< "End\n" << endl;
return 0;
}

View File

@ -34,6 +34,9 @@ Usage
\param -test \n
Suppress writing the updated files with split cyclics
\param -enableFunctionEntries \n
By default all dictionary preprocessing of fields is disabled
\*---------------------------------------------------------------------------*/
#include "argList.H"
@ -391,7 +394,12 @@ int main(int argc, char *argv[])
{
timeSelector::addOptions();
argList::addBoolOption("test");
argList::addBoolOption("test", "test only; do not change any files");
argList::addBoolOption
(
"enableFunctionEntries",
"enable expansion of dictionary directives - #include, #codeStream etc"
);
# include "addRegionOption.H"
# include "setRootCase.H"
@ -404,6 +412,7 @@ int main(int argc, char *argv[])
{
Info<< "-test option: no changes made" << nl << endl;
}
const bool enableEntries = args.optionFound("enableFunctionEntries");
Foam::word regionName = polyMesh::defaultRegion;
@ -482,6 +491,13 @@ int main(int argc, char *argv[])
IOobjectList objects(runTime, runTime.timeName());
int oldFlag = entry::disableFunctionEntries;
if (!enableEntries)
{
// By default disable dictionary expansion for fields
entry::disableFunctionEntries = 1;
}
// volFields
// ~~~~~~~~~
@ -615,6 +631,8 @@ int main(int argc, char *argv[])
thisNames,
nbrNames
);
entry::disableFunctionEntries = oldFlag;
}
return 0;

View File

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

View File

@ -0,0 +1,15 @@
EXE_INC = \
-I$(LIB_SRC)/triSurface/lnInclude \
-I$(LIB_SRC)/parallel/distributed/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/OpenFOAM/lnInclude
EXE_LIBS = \
-lmeshTools \
-lfiniteVolume \
-lOpenFOAM \
-lmeshTools \
-ltriSurface \
-ldistributed \
-lradiationModels

View File

@ -0,0 +1,54 @@
Random rndGen(653213);
// Determine mesh bounding boxes:
List<treeBoundBox> meshBb
(
1,
treeBoundBox
(
boundBox(coarseMesh.points(), false)
).extend(rndGen, 1E-3)
);
// Dummy bounds dictionary
dictionary dict;
dict.add("bounds", meshBb);
dict.add
(
"distributionType",
distributedTriSurfaceMesh::distributionTypeNames_
[
distributedTriSurfaceMesh::FROZEN
]
);
dict.add("mergeDistance", SMALL);
labelHashSet includePatches;
forAll(patches, patchI) //
{
if (!isA<processorPolyPatch>(patches[patchI]))
{
includePatches.insert(patchI);
}
}
distributedTriSurfaceMesh surfacesMesh
(
IOobject
(
"wallSurface.stl",
runTime.constant(), // directory
"triSurface", // instance
runTime, // registry
IOobject::NO_READ,
IOobject::NO_WRITE
),
triSurfaceTools::triangulate
(
patches,
includePatches
),
dict
);
//surfacesMesh.searchableSurface::write();

View File

@ -0,0 +1,97 @@
// All rays expressed as start face (local) index end end face (global)
// Pre-size by assuming a certain percentage is visible.
// Maximum lenght for dynamicList
const label maxDynListLenght = 10000;
//label lenghtCount = 0;
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
// Shoot rays from me to procI. Note that even if processor has
// 0 faces we still need to call findLine to keep calls synced.
DynamicField<point> start(coarseMesh.nFaces());
DynamicField<point> end(start.size());
DynamicList<label> startIndex(start.size());
DynamicList<label> endIndex(start.size());
const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
const pointField& remoteArea = remoteCoarseSf[procI];
const pointField& remoteFc = remoteCoarseCf[procI];
if (myFc.size()*remoteFc.size() > 0)
{
forAll(myFc, i)
{
const point& fc = myFc[i];
const vector& fA = myArea[i];
forAll(remoteFc, j)
{
if (procI != Pstream::myProcNo() || i != j)
{
const point& remFc = remoteFc[j];
const vector& remA = remoteArea[j];
const vector& d = remFc-fc;
if (((d & fA) < 0.) && ((d & remA) > 0))
{
//lenghtCount ++;
start.append(fc + 0.0001*d);
startIndex.append(i);
end.append(fc + 0.9999*d);
label globalI = globalNumbering.toGlobal(procI, j);
endIndex.append(globalI);
if (startIndex.size() > maxDynListLenght)
{
List<pointIndexHit> hitInfo(startIndex.size());
surfacesMesh.findLine
(
start,
end,
hitInfo
);
surfacesMesh.findLine(start, end, hitInfo);
forAll (hitInfo, rayI)
{
if (!hitInfo[rayI].hit())
{
rayStartFace.append(startIndex[rayI]);
rayEndFace.append(endIndex[rayI]);
}
}
//lenghtCount = 0;
start.clear();
startIndex.clear();
end.clear();
endIndex.clear();
}
}
}
}
}
}
if (!start.empty())
{
List<pointIndexHit> hitInfo(startIndex.size());
surfacesMesh.findLine
(
start,
end,
hitInfo
);
surfacesMesh.findLine(start, end, hitInfo);
forAll (hitInfo, rayI)
{
if (!hitInfo[rayI].hit())
{
rayStartFace.append(startIndex[rayI]);
rayEndFace.append(endIndex[rayI]);
}
}
}
}

View File

@ -0,0 +1,861 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ 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
viewFactorGenerator
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 "fvMesh.H"
#include "Time.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "distributedTriSurfaceMesh.H"
#include "triSurfaceTools.H"
#include "mapDistribute.H"
#include "OFstream.H"
#include "meshTools.H"
#include "plane.H"
#include "uindirectPrimitivePatch.H"
#include "DynamicField.H"
#include "IFstream.H"
#include "unitConversion.H"
#include "mathematicalConstants.H"
#include "scalarMatrices.H"
#include "CompactListList.H"
#include "labelIOList.H"
#include "labelListIOList.H"
#include "scalarListIOList.H"
#include "singleCellFvMesh.H"
#include "IOdictionary.H"
#include "fixedValueFvPatchFields.H"
using namespace Foam;
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);
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)
);
}
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];
}
}
}
// Main program:
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);
// Debug
// ~~~~~
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)
{
Info << "\nCreating single cell mesh..." << endl;
}
singleCellFvMesh coarseMesh
(
IOobject
(
mesh.name(),
runTime.timeName(),
runTime,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
finalAgglom
);
// - 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);
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];
forAll(coarseCf, faceI)
{
point cf = coarseCf[faceI];
const label coarseFaceI = coarsePatchFace[faceI];
const labelList& fineFaces = coarseToFine[coarseFaceI];
// 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()
).assign(upp.faceCentres());
SubList<point>
(
availablePoints,
upp.localPoints().size(),
upp.faceCentres().size()
).assign(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);
}
}
// - Collect remote Cf and Sf on coarse mesh //.
//------------------------------------------//
List<pointField> remoteCoarseCf(Pstream::nProcs());
List<pointField> remoteCoarseSf(Pstream::nProcs());
remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
// - Collect remote Cf and Sf on fine mesh //.
//------------------------------------------//
List<pointField> remoteFineCf(Pstream::nProcs());
List<pointField> remoteFineSf(Pstream::nProcs());
remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
// Distribute local coarse Cf and Sf for shooting rays
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Pstream::gatherList(remoteCoarseCf);
Pstream::scatterList(remoteCoarseCf);
Pstream::gatherList(remoteCoarseSf);
Pstream::scatterList(remoteCoarseSf);
// 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());
globalIndex globalNumbering(nCoarseFaces);
//- 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(), pTraits<vector>::zero);
pointField compactCoarseSf(map.constructSize(), pTraits<vector>::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).assign(localCoarseSf);
SubList<point>(compactCoarseCf, nCoarseFaces).assign(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,
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 = vector::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, 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;
}
// ************************************************************************* //