91 label localTriFaceI = 0;
93 for (
const label patchI : includePatches)
108 f.triangles(
points, nTri, triFaces);
110 forAll(triFaces, triFaceI)
112 const face&
f = triFaces[triFaceI];
121 finalAgglom[patchI][patchFaceI]
122 + coarsePatches[patchI].start()
140 rawSurface.localFaces(),
141 rawSurface.localPoints()
149 for (
const label patchI : includePatches)
153 surface.patches()[newPatchI].index() = patchI;
155 surface.patches()[newPatchI].geometricType() =
patch.type();
179 const labelList visFaces = visibleFaceFaces[faceI];
180 forAll(visFaces, faceRemote)
182 label compactI = visFaces[faceRemote];
183 const point& remoteFc = compactCf[compactI];
185 meshTools::writeOBJ(str, myFc[faceI]);
187 meshTools::writeOBJ(str, remoteFc);
189 str <<
"l " << vertI-1 <<
' ' << vertI <<
nl;
196scalar calculateViewFactorFij
205 scalar rMag =
mag(r);
209 scalar dAiMag =
mag(dAi);
210 scalar dAjMag =
mag(dAj);
214 scalar cosThetaJ =
mag(nj & r)/rMag;
215 scalar cosThetaI =
mag(ni & r)/rMag;
219 (cosThetaI*cosThetaJ*dAjMag*dAiMag)
220 /(
sqr(rMag)*constant::mathematical::pi)
230void insertMatrixElements
233 const label fromProcI,
239 forAll(viewFactors, faceI)
242 const labelList& globalFaces = globalFaceFaces[faceI];
244 label globalI = globalNumbering.
toGlobal(fromProcI, faceI);
247 matrix[globalI][globalFaces[i]] = vf[i];
255int main(
int argc,
char *argv[])
259 "Calculate view factors from face agglomeration array."
260 " The finalAgglom generated by faceAgglomerate utility."
276 IOobject::MUST_READ_IF_MODIFIED,
281 const word viewFactorWall(
"viewFactorWall");
283 const bool writeViewFactors =
284 viewFactorDict.getOrDefault(
"writeViewFactorMatrix",
false);
286 const bool dumpRays =
287 viewFactorDict.getOrDefault(
"dumpRays",
false);
289 const label
debug = viewFactorDict.getOrDefault<label>(
"debug", 0);
297 mesh.facesInstance(),
310 Pout <<
"\nCreating single cell mesh..." <<
endl;
317 "coarse:" +
mesh.name(),
329 Pout <<
"\nCreated single cell mesh..." <<
endl;
336 label nCoarseFaces = 0;
337 label nFineFaces = 0;
343 for (
const label patchi : viewFactorsPatches)
345 nCoarseFaces += coarsePatches[patchi].
size();
346 nFineFaces +=
patches[patchi].size();
350 label totalNCoarseFaces = nCoarseFaces;
354 if (Pstream::master())
356 Info <<
"\nTotal number of coarse faces: "<< totalNCoarseFaces <<
endl;
359 if (Pstream::master() && debug)
361 Pout <<
"\nView factor patches included in the calculation : "
362 << viewFactorsPatches <<
endl;
373 for (
const label
patchID : viewFactorsPatches)
380 if (agglom.
size() > 0)
382 label nAgglom =
max(agglom)+1;
385 coarseMesh.patchFaceMap()[
patchID];
388 coarseMesh.Cf().boundaryField()[
patchID];
390 coarseMesh.Sf().boundaryField()[
patchID];
394 point cf = coarseCf[faceI];
396 const label coarseFaceI = coarsePatchFace[faceI];
397 const labelList& fineFaces = coarseToFine[coarseFaceI];
398 const label agglomI =
410 upp.faceCentres().size()
411 + upp.localPoints().size()
417 upp.faceCentres().size()
418 ) = upp.faceCentres();
423 upp.localPoints().size(),
424 upp.faceCentres().size()
425 ) = upp.localPoints();
429 forAll(availablePoints, iPoint)
431 point cfFine = availablePoints[iPoint];
432 if (
mag(cfFine-cfo) < dist)
434 dist =
mag(cfFine-cfo);
439 point sf = coarseSf[faceI];
441 localCoarseSf.append(sf);
442 localAgg.append(agglomI);
455 remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
456 remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
457 remoteCoarseAgg[Pstream::myProcNo()] = localAgg;
459 Pstream::allGatherList(remoteCoarseCf);
460 Pstream::allGatherList(remoteCoarseSf);
461 Pstream::allGatherList(remoteCoarseAgg);
483 labelList nVisibleFaceFaces(nCoarseFaces, Zero);
487 nVisibleFaceFaces[rayStartFace[i]]++;
492 label nViewFactors = 0;
493 forAll(nVisibleFaceFaces, faceI)
495 visibleFaceFaces[faceI].
setSize(nVisibleFaceFaces[faceI]);
496 nViewFactors += nVisibleFaceFaces[faceI];
513 nVisibleFaceFaces = 0;
516 label faceI = rayStartFace[i];
517 label compactI = rayEndFace[i];
518 visibleFaceFaces[faceI][nVisibleFaceFaces[faceI]++] = compactI;
526 pointField compactCoarseCf(map.constructSize(), Zero);
527 pointField compactCoarseSf(map.constructSize(), Zero);
539 forAll(viewFactorsPatches, i)
541 label
patchID = viewFactorsPatches[i];
544 if (agglom.
size() > 0)
546 label nAgglom =
max(agglom)+1;
549 coarseMesh.patchFaceMap()[
patchID];
551 forAll(coarseToFine, coarseI)
553 compactPatchId.append(
patchID);
557 const label coarseFaceI = coarsePatchFace[coarseI];
558 const labelList& fineFaces = coarseToFine[coarseFaceI];
566 coarseToFine[coarseFaceI]
571 coarseToFine[coarseFaceI]
578 map.distribute(compactCoarseSf);
579 map.distribute(compactCoarseCf);
580 map.distribute(compactFineCf);
581 map.distribute(compactFineSf);
583 map.distribute(compactPatchId);
591 runTime.path()/
"allVisibleFaces.obj",
593 remoteCoarseCf[Pstream::myProcNo()],
607 mesh.facesInstance(),
616 label totalPatches = coarsePatches.
size();
628 if (Pstream::master())
630 Info<<
"\nCalculating view factors..." <<
endl;
633 if (
mesh.nSolutionD() == 3)
635 forAll(localCoarseSf, coarseFaceI)
637 const List<point>& localFineSf = compactFineSf[coarseFaceI];
639 const List<point>& localFineCf = compactFineCf[coarseFaceI];
640 const label fromPatchId = compactPatchId[coarseFaceI];
641 patchArea[fromPatchId] +=
mag(Ai);
643 const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
645 forAll(visCoarseFaces, visCoarseFaceI)
647 F[coarseFaceI].setSize(visCoarseFaces.
size());
648 label compactJ = visCoarseFaces[visCoarseFaceI];
649 const List<point>& remoteFineSj = compactFineSf[compactJ];
650 const List<point>& remoteFineCj = compactFineCf[compactJ];
652 const label toPatchId = compactPatchId[compactJ];
657 const vector& dAi = localFineSf[i];
658 const vector& dCi = localFineCf[i];
662 const vector& dAj = remoteFineSj[j];
663 const vector& dCj = remoteFineCj[j];
665 scalar dIntFij = calculateViewFactorFij
676 F[coarseFaceI][visCoarseFaceI] = Fij/
mag(Ai);
677 sumViewFactorPatch[fromPatchId][toPatchId] += Fij;
681 else if (
mesh.nSolutionD() == 2)
694 scalar wideBy2 = (box.
span() & emptyDir)*2.0;
696 forAll(localCoarseSf, coarseFaceI)
698 const vector& Ai = localCoarseSf[coarseFaceI];
699 const vector& Ci = localCoarseCf[coarseFaceI];
701 vector R1i = Ci + (
mag(Ai)/wideBy2)*(Ain ^ emptyDir);
702 vector R2i = Ci - (
mag(Ai)/wideBy2)*(Ain ^ emptyDir) ;
704 const label fromPatchId = compactPatchId[coarseFaceI];
705 patchArea[fromPatchId] +=
mag(Ai);
707 const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
708 forAll(visCoarseFaces, visCoarseFaceI)
710 F[coarseFaceI].setSize(visCoarseFaces.
size());
711 label compactJ = visCoarseFaces[visCoarseFaceI];
712 const vector& Aj = compactCoarseSf[compactJ];
713 const vector& Cj = compactCoarseCf[compactJ];
715 const label toPatchId = compactPatchId[compactJ];
718 vector R1j = Cj + (
mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
719 vector R2j = Cj - (
mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
721 scalar d1 =
mag(R1i - R2j);
722 scalar d2 =
mag(R2i - R1j);
723 scalar s1 =
mag(R1i - R1j);
724 scalar s2 =
mag(R2i - R2j);
726 scalar Fij =
mag((d1 + d2) - (s1 + s2))/(4.0*
mag(Ai)/wideBy2);
728 F[coarseFaceI][visCoarseFaceI] = Fij;
729 sumViewFactorPatch[fromPatchId][toPatchId] += Fij*
mag(Ai);
734 if (Pstream::master())
736 Info <<
"Writing view factor matrix..." <<
endl;
746 if (Pstream::master() && debug)
748 forAll(viewFactorsPatches, i)
750 label patchI = viewFactorsPatches[i];
751 forAll(viewFactorsPatches, i)
753 label patchJ = viewFactorsPatches[i];
754 Info <<
"F" << patchI << patchJ <<
": "
755 << sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
762 if (writeViewFactors)
769 mesh.time().timeName(),
781 forAll(viewFactorsPatches, i)
783 label
patchID = viewFactorsPatches[i];
785 if (agglom.
size() > 0)
787 label nAgglom =
max(agglom)+1;
790 coarseMesh.patchFaceMap()[
patchID];
792 forAll(coarseToFine, coarseI)
794 const scalar Fij =
sum(
F[compactI]);
795 const label coarseFaceID = coarsePatchFace[coarseI];
796 const labelList& fineFaces = coarseToFine[coarseFaceID];
799 const label faceID = fineFaces[fineId];
806 viewFactorField.write();
813 labelList compactToGlobal(map.constructSize());
816 for (label i = 0; i < globalNumbering.
localSize(); i++)
818 compactToGlobal[i] = globalNumbering.
toGlobal(i);
824 const Map<label>& localToCompactMap = compactMap[procI];
828 compactToGlobal[*iter] = globalNumbering.
toGlobal
841 forAll(globalFaceFaces, faceI)
846 visibleFaceFaces[faceI]
855 mesh.facesInstance(),
861 std::move(globalFaceFaces)
863 IOglobalFaceFaces.write();
871 mesh.facesInstance(),
reduce(hasMovingMesh, orOp< bool >())
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
IOmapDistribute is derived from mapDistribute and IOobject to give the mapDistribute automatic IO fun...
Defines the attributes of an object for which implicit objectRegistry management is supported,...
void setSize(const label n)
Alias for resize()
void append(const T &val)
Append an element at the end of the list.
A HashTable to objects of type <T> with a label key.
Output to file stream, using an OSstream.
virtual const fileName & name() const
Get the name of the stream.
A list of faces which address into the list of points.
const Field< point_type > & points() const noexcept
Return reference to global points.
A List obtained as a section of another List.
A List with indirect addressing. Like IndirectList but does not store addressing.
void size(const label n)
Older name for setAddressableSize.
label size() const noexcept
The number of elements in the list.
A bounding box defined in terms of min/max extrema points.
vector span() const
The bounding box span (from minimum to maximum)
A face is a list of labels corresponding to mesh vertices.
A class for handling file names.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
label localSize() const
My local size.
label toGlobal(const label i) const
From local to global index.
A triFace with additional (region) index.
Class containing processor-to-processor mapping information.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
label start() const noexcept
The start label of boundary faces in the polyMesh face list.
Mesh consisting of general polyhedral cells.
A patch is a list of labels that address the faces in the global face list.
fvMesh as subset of other mesh. Consists of one cell and all original bounday faces....
Triangulated surface description with patch information.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
const polyBoundaryMesh & patches
volVectorField F(fluid.F())
int debug
Static debugging option.
const wordList surface
Standard surface field types (scalar, vector, tensor, etc)
const std::string patch
OpenFOAM patch number as a std::string.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
IntListType renumber(const labelUList &oldToNew, const IntListType &input)
Renumber the values (not the indices) of a list.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
constexpr char nl
The newline '\n' character (0x0a)
labelList triSurfaceToAgglom(5 *nFineFaces)
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Unit conversion functions.