45 void Foam::faceReflecting::initialise(
const dictionary& coeffs)
48 forAll(qreflective_, bandI)
58 mesh_.time().timeName(),
70 if (mesh_.nSolutionD() == 3)
72 nRay_ = 4*nPhi_*nTheta_;
73 refDiscAngles_.resize(nRay_);
74 const scalar deltaPhi =
pi/(2.0*nPhi_);
75 const scalar deltaTheta =
pi/nTheta_;
77 for (label
n = 1;
n <= nTheta_;
n++)
79 for (label m = 1; m <= 4*nPhi_; m++)
81 const scalar thetai = (2*
n - 1)*deltaTheta/2.0;
82 const scalar phii = (2*m - 1)*deltaPhi/2.0;
88 refDiscAngles_[rayI++] =
95 else if (mesh_.nSolutionD() == 2)
98 refDiscAngles_.resize(nRay_);
101 const scalar deltaPhi =
pi/(2.0*nPhi_);
102 for (label m = 1; m <= 4*nPhi_; m++)
104 const scalar phii = (2*m - 1)*deltaPhi/2.0;
111 refDiscAngles_[rayI++] =
118 <<
"The reflected rays are available in 2D or 3D " 122 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
124 const radiation::boundaryRadiationProperties& boundaryRadiation =
128 globalIndex globalNumbering(mesh_.nFaces());
133 DynamicList<point> dynCf;
134 DynamicList<vector> dynNf;
135 DynamicList<label> dynFacesI;
138 const polyPatch& pp =
patches[patchI];
141 if (!pp.coupled() && !isA<cyclicAMIPolyPatch>(pp))
143 const tmp<scalarField> tt =
144 boundaryRadiation.transmissivity(patchI);
146 const tmp<scalarField>
tr =
147 boundaryRadiation.specReflectivity(patchI);
149 const tmp<scalarField> ta =
150 boundaryRadiation.absorptivity(patchI);
169 dynFacesI.append(faceI + pp.start());
170 dynCf.append(cf[faceI]);
171 dynNf.append(
n[faceI]);
178 (r[faceI] > 0 && t[faceI] == 0) ||
179 (t[faceI] == 0 && a[faceI] > 0 && r[faceI] == 0)
182 includePatches_.insert(patchI);
188 shootFacesIds_.reset(
new labelList(dynFacesI));
200 treeBoundBox(mesh_.points()).extend(
rndGen, 1
e-3)
214 dict.add(
"mergeDistance", SMALL);
219 mesh_.boundaryMesh(),
226 new distributedTriSurfaceMesh
230 "reflectiveSurface.stl",
231 mesh_.time().constant(),
244 surfacesMesh_->searchableSurface::write();
249 void Foam::faceReflecting::calculate()
251 const radiation::boundaryRadiationProperties& boundaryRadiation =
256 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
258 const fvBoundaryMesh& fvPatches = mesh_.boundary();
260 label nBands = spectralDistribution_.size();
264 const vector sunDir = directHitFaces_.direction();
265 const labelList& directHits = directHitFaces_.rayStartFaces();
267 globalIndex globalNumbering(mesh_.nFaces());
269 Map<label> refFacesDirIndex;
274 const polyPatch& pp =
patches[patchI];
276 if (!pp.coupled() && !isA<cyclicAMIPolyPatch>(pp))
278 const tmp<scalarField>
tr =
279 boundaryRadiation.specReflectivity(patchI);
286 label globalID = faceI + pp.start();
288 if (r[faceI] > 0.0 && directHits.found(globalID))
291 sunDir + 2.0*(-sunDir &
n[faceI]) *
n[faceI];
296 forAll(refDiscAngles_, iDisc)
298 scalar dotProd = refDir & refDiscAngles_[iDisc];
308 if (refDisDirsIndex[rayIndex] == -1)
310 refDisDirsIndex[rayIndex] = 1;
314 refFacesDirIndex.insert
316 globalNumbering.toGlobal(globalID),
332 const scalar maxBounding =
333 returnReduce(5.0*mesh_.bounds().mag(), maxOp<scalar>());
339 DynamicField<point> start(nFaces);
340 DynamicField<point>
end(start.size());
341 DynamicList<label> startIndex(start.size());
342 DynamicField<label> dirStartIndex(start.size());
347 for (; i < Cfs_->size(); i++)
349 const point& fc = Cfs_()[i];
351 const vector nf = Nfs_()[i];
353 const label myFaceId = shootFacesIds_()[i];
355 forAll(refDisDirsIndex, dirIndex)
357 if (refDisDirsIndex[dirIndex] > -1)
359 if ((nf & refDiscAngles_[dirIndex]) > 0)
365 startIndex.append(myFaceId);
366 dirStartIndex.append(dirIndex);
376 List<pointIndexHit> hitInfo(startIndex.size());
378 surfacesMesh_->findLine(start,
end, hitInfo);
382 autoPtr<mapDistribute> mapPtr
384 surfacesMesh_->localQueries
390 const mapDistribute& map = mapPtr();
392 PtrList<List<scalarField>> patchr(
patches.size());
393 PtrList<List<scalarField>> patcha(
patches.size());
399 new List<scalarField>(nBands)
405 new List<scalarField>(nBands)
412 const polyPatch& pp =
patches[patchi];
414 if (!pp.coupled() && !isA<cyclicAMIPolyPatch>(pp))
416 for (label bandI = 0; bandI < nBands; bandI++)
418 patchr[patchi][bandI] =
419 boundaryRadiation.specReflectivity
426 patcha[patchi][bandI] =
427 boundaryRadiation.absorptivity
437 List<scalarField> r(nBands);
438 for (label bandI = 0; bandI < nBands; bandI++)
440 r[bandI].setSize(triangleIndex.size());
442 labelList refDirIndex(triangleIndex.size(), -1);
443 labelList refIndex(triangleIndex.size(), -1);
447 label trii = triangleIndex[i];
448 label facei = mapTriToGlobal_[trii];
449 label patchI =
patches.whichPatch(facei);
450 const polyPatch& pp =
patches[patchI];
451 label localFaceI = pp.whichFace(facei);
454 if (refFacesDirIndex.found(globalFace))
456 refDirIndex[i] = refFacesDirIndex.find(globalFace)();
457 refIndex[i] = globalFace;
459 for (label bandI = 0; bandI < nBands; bandI++)
461 r[bandI][i] = patchr[patchI][bandI][localFaceI];
464 map.reverseDistribute(hitInfo.size(), refDirIndex);
465 map.reverseDistribute(hitInfo.size(), refIndex);
466 for (label bandI = 0; bandI < nBands; bandI++)
468 map.reverseDistribute(hitInfo.size(), r[bandI]);
471 for (label bandI = 0; bandI < nBands; bandI++)
474 qreflective_[bandI].boundaryFieldRef();
478 const vector qPrim(solarCalc_.directSolarRad()*solarCalc_.direction());
484 if (hitInfo[rayI].hit())
488 dirStartIndex[rayI] == refDirIndex[rayI]
489 && refFacesDirIndex.found(refIndex[rayI])
492 for (label bandI = 0; bandI < nBands; bandI++)
495 qreflective_[bandI].boundaryFieldRef();
497 label startFaceId = startIndex[rayI];
498 label startPatchI =
patches.whichPatch(startFaceId);
500 const polyPatch& ppStart =
patches[startPatchI];
501 label localStartFaceI = ppStart.whichFace(startFaceId);
503 scalar a = patcha[startPatchI][bandI][localStartFaceI];
507 vector rayIn = refDiscAngles_[dirStartIndex[rayI]];
511 qrefBf[startPatchI][localStartFaceI] +=
516 *spectralDistribution_[bandI]
520 & nStart[localStartFaceI]
530 dirStartIndex.clear();
536 Foam::faceReflecting::faceReflecting
546 nTheta_(
dict.subDict(
"reflecting").getOrDefault<label>(
"nTheta", 10)),
547 nPhi_(
dict.subDict(
"reflecting").getOrDefault<label>(
"nPhi", 10)),
550 spectralDistribution_(spectralDistribution),
551 qreflective_(spectralDistribution_.size()),
552 directHitFaces_(directHiyFaces),
static const Enum< distributionType > distributionTypeNames_
const triSurface localSurface
A solar calculator model providing models for the solar direction and solar loads.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
static const boundaryRadiationProperties & New(const fvMesh &mesh, Args &&... args)
Get existing or create a new MeshObject.
Ignore writing from objectRegistry::writeObject()
static int myProcNo(const label communicator=worldComm)
Number of this process (starting from masterNo() = 0)
SubField is a Field obtained as a section of another Field, without its own allocation. SubField is derived from a SubList rather than a List.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
#define forAll(list, i)
Loop across all elements in list.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
vectorField pointField
pointField is a vectorField.
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
dimensionedScalar cos(const dimensionedScalar &ds)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
constexpr scalar pi(M_PI)
dimensionedScalar tr(const dimensionedSphericalTensor &dt)
constexpr scalar piByTwo(0.5 *M_PI)
errorManip< error > abort(error &err)
Helper class to calculate visible faces for global, sun-like illumination.
dimensionedScalar sin(const dimensionedScalar &ds)
int debug
Static debugging option.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
defineTypeNameAndDebug(combustionModel, 0)
void correct()
Correct reflected flux.
dimensionedScalar pow3(const dimensionedScalar &ds)
vector point
Point is a vector.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Mesh data needed to do the Finite Volume discretisation.
const polyBoundaryMesh & patches
Automatically write from objectRegistry::writeObject()
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
void clear()
Clear the PtrList. Delete allocated entries and set size to zero.
Field< vector > vectorField
Specialisation of Field<T> for vector.
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
List< label > labelList
A List of labels.
static void mapCombineReduce(Container &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) applying cop to inplace combine map values from different processo...
bool returnReduceOr(const bool value, const label comm=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
List< treeBoundBox > meshBb(1, treeBoundBox(coarseMesh.points()).extend(rndGen, 1e-3))
static void listCombineReduce(List< T > &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
After completion all processors have the same data.
static constexpr const zero Zero
Global zero (0)