97 #include "CGAL/version.h" 98 #if defined(CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1050211000) 99 #define BOOST_BIND_GLOBAL_PLACEHOLDERS 101 #pragma clang diagnostic ignored "-Wbitwise-instead-of-logical" 104 #include <CGAL/Simple_cartesian.h> 105 #include <CGAL/AABB_tree.h> 106 #include <CGAL/AABB_traits.h> 107 #include <CGAL/AABB_triangle_primitive.h> 108 #include <CGAL/Surface_mesh.h> 110 typedef CGAL::Simple_cartesian<double>
K;
111 typedef K::Point_3
Point;
112 typedef K::Direction_3 Vector3C;
116 typedef std::vector<Triangle>::iterator Iterator;
117 typedef CGAL::AABB_triangle_primitive<K, Iterator> Primitive;
118 typedef CGAL::AABB_traits<K, Primitive> AABB_triangle_traits;
119 typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
120 typedef boost::optional
122 Tree::Intersection_and_primitive_id<Segment>::Type
123 > Segment_intersection;
127 using namespace Foam;
147 label localTriFaceI = 0;
149 for (
const label patchI : includePatches)
164 f.triangles(
points, nTri, triFaces);
166 forAll(triFaces, triFaceI)
168 const face&
f = triFaces[triFaceI];
177 finalAgglom[patchI][patchFaceI]
178 + coarsePatches[patchI].start()
212 for (
const label proci : globalPointIdx.allProcs())
214 const label offset = globalPointIdx.localStart(proci);
221 : globalSurfaceTris.slice(globalFaceIdx.range(proci))
234 std::move(globalSurfaceTris),
235 std::move(globalSurfacePoints)
246 for (
const label patchI : includePatches)
250 surface.patches()[newPatchI].index() = patchI;
252 surface.patches()[newPatchI].geometricType() =
patch.type();
276 const labelList visFaces = visibleFaceFaces[faceI];
277 forAll(visFaces, faceRemote)
279 label compactI = visFaces[faceRemote];
280 const point& remoteFc = compactCf[compactI];
286 str <<
"l " << vertI-1 <<
' ' << vertI <<
nl;
293 scalar calculateViewFactorFij2AI
302 scalar rMag =
mag(r);
306 scalar dAiMag =
mag(dAi);
307 scalar dAjMag =
mag(dAj);
311 scalar cosThetaJ =
mag(nj & r)/rMag;
312 scalar cosThetaI =
mag(ni & r)/rMag;
316 (cosThetaI*cosThetaJ*dAjMag*dAiMag)
327 void insertMatrixElements
330 const label fromProcI,
336 forAll(viewFactors, faceI)
339 const labelList& globalFaces = globalFaceFaces[faceI];
341 label globalI = globalNumbering.
toGlobal(fromProcI, faceI);
344 matrix[globalI][globalFaces[i]] = vf[i];
372 r = (
alpha*magSi)*di;
381 pi[i] = ci +
p[i]*(magSi/2)*di;
387 pj[i] = cj +
p[i]*(magSj/2)*dj;
398 r = (
alpha*magSi)*di;
411 dIntFij *= (magSi/2) * (magSj/2);
420 int main(
int argc,
char *argv[])
424 "Calculate view factors from face agglomeration array." 425 " The finalAgglom generated by faceAgglomerate utility." 446 const word viewFactorWall(
"viewFactorWall");
448 const bool writeViewFactors =
449 viewFactorDict.getOrDefault(
"writeViewFactorMatrix",
false);
451 const bool dumpRays =
452 viewFactorDict.getOrDefault(
"dumpRays",
false);
454 const label
debug = viewFactorDict.getOrDefault<label>(
"debug", 0);
456 const scalar GaussQuadTol =
457 viewFactorDict.getOrDefault<scalar>(
"GaussQuadTol", 0.01);
459 const scalar distTol =
460 viewFactorDict.getOrDefault<scalar>(
"distTol", 8);
465 const scalar intTol =
466 viewFactorDict.getOrDefault<scalar>(
"intTol", 1
e-2);
468 bool useAgglomeration(
true);
490 for (label patchi=0; patchi <
patches.
size(); patchi++)
494 useAgglomeration =
false;
502 Pout <<
"\nCreating single cell mesh..." <<
endl;
521 Pout <<
"\nCreated single cell mesh..." <<
endl;
528 label nCoarseFaces = 0;
529 label nFineFaces = 0;
533 for (
const label patchi : viewFactorsPatches)
535 nCoarseFaces += coarsePatches[patchi].
size();
540 Info<<
"\nTotal number of coarse faces: " 546 Pout <<
"\nView factor patches included in the calculation : " 547 << viewFactorsPatches <<
endl;
558 for (
const label
patchID : viewFactorsPatches)
565 if (agglom.
size() > 0)
567 label nAgglom =
max(agglom)+1;
570 coarseMesh.patchFaceMap()[
patchID];
573 coarseMesh.Cf().boundaryField()[
patchID];
575 coarseMesh.Sf().boundaryField()[
patchID];
579 point cf = coarseCf[faceI];
581 const label coarseFaceI = coarsePatchFace[faceI];
582 const labelList& fineFaces = coarseToFine[coarseFaceI];
583 const label agglomI =
595 upp.faceCentres().size()
596 + upp.localPoints().size()
602 upp.faceCentres().
size()
603 ) = upp.faceCentres();
608 upp.localPoints().size(),
609 upp.faceCentres().size()
610 ) = upp.localPoints();
614 forAll(availablePoints, iPoint)
616 point cfFine = availablePoints[iPoint];
617 if (
mag(cfFine-cfo) < dist)
619 dist =
mag(cfFine-cfo);
624 point sf = coarseSf[faceI];
626 localCoarseSf.append(sf);
627 localAgg.append(agglomI);
682 nVisibleFaceFaces[rayStartFace[i]]++;
687 label nViewFactors = 0;
688 forAll(nVisibleFaceFaces, faceI)
690 visibleFaceFaces[faceI].
setSize(nVisibleFaceFaces[faceI]);
691 nViewFactors += nVisibleFaceFaces[faceI];
708 nVisibleFaceFaces = 0;
711 label faceI = rayStartFace[i];
712 label compactI = rayEndFace[i];
713 visibleFaceFaces[faceI][nVisibleFaceFaces[faceI]++] = compactI;
738 forAll(viewFactorsPatches, i)
740 label
patchID = viewFactorsPatches[i];
743 if (agglom.
size() > 0)
745 label nAgglom =
max(agglom)+1;
748 coarseMesh.patchFaceMap()[
patchID];
752 forAll(coarseToFine, coarseI)
754 compactPatchId.append(
patchID);
758 label startFace = pp.
start();
763 faces[coarseI + startFace]
766 const label coarseFaceI = coarsePatchFace[coarseI];
767 const labelList& fineFaces = coarseToFine[coarseFaceI];
772 compactPoints.append(locPoints);
777 coarseToFine[coarseFaceI]
782 coarseToFine[coarseFaceI]
796 map.distribute(compactCoarseSf);
797 map.distribute(compactCoarseCf);
798 map.distribute(compactFineCf);
799 map.distribute(compactFineSf);
800 map.distribute(compactPoints);
801 map.distribute(compactPatchId);
832 const label totalPatches =
842 Info<<
"\nCalculating view factors..." <<
endl;
854 GaussPoints[2][0] = 0;
860 GaussPoints[3][1] = -GaussPoints[3][0];
862 GaussPoints[3][3] = -GaussPoints[3][2];
865 GaussPoints[4][0] = 0;
867 GaussPoints[4][2] = -GaussPoints[4][1];
869 GaussPoints[4][4] = -GaussPoints[4][3];
877 GaussWeights[1][0] = 1;
878 GaussWeights[1][1] = 1;
881 GaussWeights[2][0] = 8.0/9.0;
882 GaussWeights[2][1] = 5.0/9.0;
883 GaussWeights[2][2] = 5.0/9.0;
886 GaussWeights[3][0] = (18.0 +
std::sqrt(30))/36.0;
887 GaussWeights[3][1] = (18.0 +
std::sqrt(30))/36.0;
888 GaussWeights[3][2] = (18.0 -
std::sqrt(30))/36.0;
889 GaussWeights[3][3] = (18.0 -
std::sqrt(30))/36.0;
892 GaussWeights[4][0] = 128.0/225.0;
893 GaussWeights[4][1] = (322.0 + 13.0*
std::sqrt(70))/900.0;
894 GaussWeights[4][2] = (322.0 + 13.0*
std::sqrt(70))/900.0;
895 GaussWeights[4][3] = (322.0 - 13.0*
std::sqrt(70))/900.0;
896 GaussWeights[4][4] = (322.0 - 13.0*
std::sqrt(70))/900.0;
898 const label maxQuadOrder = 5;
902 forAll(localCoarseSf, coarseFaceI)
904 const List<point>& localFineSf = compactFineSf[coarseFaceI];
906 const List<point>& localFineCf = compactFineCf[coarseFaceI];
907 const label fromPatchId = compactPatchId[coarseFaceI];
909 const List<point>& lPoints = compactPoints[coarseFaceI];
911 patchArea[fromPatchId] +=
mag(Ai);
913 const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
915 forAll(visCoarseFaces, visCoarseFaceI)
918 F2LI[coarseFaceI].setSize(visCoarseFaces.
size());
919 label compactJ = visCoarseFaces[visCoarseFaceI];
920 const List<point>& remoteFineSj = compactFineSf[compactJ];
921 const List<point>& remoteFineCj = compactFineCf[compactJ];
923 const List<point>& rPointsCj = compactPoints[compactJ];
925 const label toPatchId = compactPatchId[compactJ];
936 const vector& dCi = localFineCf[i];
945 const vector& dCj = remoteFineCj[j];
947 const scalar dist =
mag(dCi - dCj)/((dAi + dAj)/2);
963 const vector& dAi = localFineSf[i];
964 const vector& dCi = localFineCf[i];
968 const vector& dAj = remoteFineSj[j];
969 const vector& dCj = remoteFineCj[j];
971 scalar dIntFij = calculateViewFactorFij2AI
982 F2LI[coarseFaceI][visCoarseFaceI] = F2AIij/
mag(Ai);
987 label nLocal = lPoints.
size();
988 label nRemote = rPointsCj.
size();
991 scalar oldEToeInt = 0;
992 for (label gi=0; gi < maxQuadOrder; gi++)
995 for(label i=0; i<nLocal; i++)
1005 si = lPoints[i] - lPoints[nLocal-1];
1006 ci = (lPoints[i] + lPoints[nLocal-1])/2;
1010 si = lPoints[i] - lPoints[i-1];
1011 ci = (lPoints[i] + lPoints[i-1])/2;
1014 for(label j=0; j<nRemote; j++)
1018 sj = rPointsCj[j]-rPointsCj[nRemote-1];
1019 cj = (rPointsCj[j]+rPointsCj[nRemote-1])/2;
1023 sj = rPointsCj[j] - rPointsCj[j-1];
1024 cj = (rPointsCj[j] + rPointsCj[j-1])/2;
1028 scalar magSi =
mag(si);
1029 scalar magSj =
mag(sj);
1030 scalar cosij = (si & sj)/(magSi * magSj);
1035 label quadOrder = gi;
1065 mag(F2LIij) > ROOTVSMALL
1066 ? (F2LIij-oldEToeInt)/F2LIij
1072 (
mag(err) < GaussQuadTol && gi > 0)
1073 || gi == maxQuadOrder-1
1076 F2LI[coarseFaceI][visCoarseFaceI] =
1082 oldEToeInt = F2LIij;
1087 sumViewFactorPatch[fromPatchId][toPatchId] +=
1088 F2LI[coarseFaceI][visCoarseFaceI]*
mag(Ai);
1105 scalar wideBy2 = (box.
span() & emptyDir)*2.0;
1107 forAll(localCoarseSf, coarseFaceI)
1109 const vector& Ai = localCoarseSf[coarseFaceI];
1110 const vector& Ci = localCoarseCf[coarseFaceI];
1112 vector R1i = Ci + (
mag(Ai)/wideBy2)*(Ain ^ emptyDir);
1113 vector R2i = Ci - (
mag(Ai)/wideBy2)*(Ain ^ emptyDir) ;
1115 const label fromPatchId = compactPatchId[coarseFaceI];
1116 patchArea[fromPatchId] +=
mag(Ai);
1118 const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
1119 forAll(visCoarseFaces, visCoarseFaceI)
1121 F2LI[coarseFaceI].setSize(visCoarseFaces.
size());
1122 label compactJ = visCoarseFaces[visCoarseFaceI];
1123 const vector& Aj = compactCoarseSf[compactJ];
1124 const vector& Cj = compactCoarseCf[compactJ];
1126 const label toPatchId = compactPatchId[compactJ];
1129 vector R1j = Cj + (
mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
1130 vector R2j = Cj - (
mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
1132 scalar d1 =
mag(R1i - R2j);
1133 scalar d2 =
mag(R2i - R1j);
1134 scalar s1 =
mag(R1i - R1j);
1135 scalar s2 =
mag(R2i - R2j);
1137 scalar Fij =
mag((d1 + d2) - (s1 + s2))/(4.0*
mag(Ai)/wideBy2);
1139 F2LI[coarseFaceI][visCoarseFaceI] = Fij;
1140 sumViewFactorPatch[fromPatchId][toPatchId] += Fij*
mag(Ai);
1147 <<
" View factors are not available in 1D " 1159 forAll(viewFactorsPatches, i)
1161 label patchI = viewFactorsPatches[i];
1162 for (label j=i; j<viewFactorsPatches.size(); j++)
1164 label patchJ = viewFactorsPatches[j];
1166 Info <<
"F" << patchI << patchJ <<
": " 1167 << sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
1173 if (writeViewFactors)
1177 Info <<
"Writing view factor matrix..." <<
endl;
1197 forAll(viewFactorsPatches, i)
1199 label
patchID = viewFactorsPatches[i];
1201 if (agglom.
size() > 0)
1203 label nAgglom =
max(agglom)+1;
1206 coarseMesh.patchFaceMap()[
patchID];
1208 forAll(coarseToFine, coarseI)
1210 const scalar FiSum =
sum(F2LI[compactI]);
1212 const label coarseFaceID = coarsePatchFace[coarseI];
1213 const labelList& fineFaces = coarseToFine[coarseFaceID];
1214 forAll(fineFaces, fineId)
1216 const label faceID = fineFaces[fineId];
1217 vfbf[
patchID][faceID] = FiSum;
1223 viewFactorField.write();
1230 labelList compactToGlobal(map.constructSize());
1233 for (label i = 0; i < globalNumbering.
localSize(); i++)
1235 compactToGlobal[i] = globalNumbering.
toGlobal(i);
1239 forAll(compactMap, procI)
1241 const Map<label>& localToCompactMap = compactMap[procI];
1245 compactToGlobal[*iter] = globalNumbering.
toGlobal 1258 forAll(globalFaceFaces, faceI)
1263 visibleFaceFaces[faceI]
1278 std::move(globalFaceFaces)
1280 IOglobalFaceFaces.
write();
Different types of constants.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
virtual const fileName & name() const
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
static void addNote(const string ¬e)
Add extra notes for the usage information.
void size(const label n)
Older name for setAddressableSize.
const surfaceVectorField & Sf() const
Return cell face area vectors.
fvMesh as subset of other mesh. Consists of one cell and all original bounday faces. Useful when manipulating boundary data. Single internal cell only needed to be able to manipulate in a standard way.
fileName path() const
Return path.
static dimensioned< Type > getOrDefault(const word &name, const dictionary &dict, const dimensionSet &dims=dimless, const Type &deflt=Type(Zero))
Construct dimensioned from dictionary, with default value.
A class for handling file names.
IOmapDistribute is derived from mapDistribute and IOobject to give the mapDistribute automatic IO fun...
dimensionedScalar log(const dimensionedScalar &ds)
errorManipArg< error, int > exit(error &err, const int errNo=1)
const fileName & facesInstance() const
Return the current instance directory for faces.
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
A 1D vector of objects of type <T> with a fixed length <N>.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const surfaceVectorField & Cf() const
Return face centres as surfaceVectorField.
void append(const T &val)
Append an element at the end of the list.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
IntListType renumber(const labelUList &oldToNew, const IntListType &input)
Renumber the values (not the indices) of a list.
Output to file stream, using an OSstream.
constexpr char nl
The newline '\n' character (0x0a)
label start() const noexcept
The start label of boundary faces in the polyMesh face list.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Dispatch tag. Construct 'one-sided' from local sizes, using gather but no broadcast.
std::vector< Triangle > triangles
A bounding box defined in terms of min/max extrema points.
CGAL::Triangle_3< K > Triangle
Ignore writing from objectRegistry::writeObject()
const dimensionSet dimless
Dimensionless.
static int myProcNo(const label communicator=worldComm)
Number of this process (starting from masterNo() = 0)
const Time & time() const
Return the top-level database.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
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.
CGAL::Exact_predicates_exact_constructions_kernel K
void setSize(const label n)
Dummy function, to make FixedList consistent with List.
static void broadcast(Type &value, const label comm=UPstream::worldComm)
Broadcast content (contiguous or non-contiguous) to all processes in communicator.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
labelList indices(const wordRe &matcher, const bool useGroups=true) const
Return (sorted) patch indices for all matches.
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
static void allGatherList(List< T > &values, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Gather data, but keep individual values separate. Uses linear/tree communication. ...
A list of faces which address into the list of points.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
A List obtained as a section of another List.
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator) is 1 for serial run.
void setSize(const label n)
Alias for resize()
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
labelList triSurfaceToAgglom(5 *nFineFaces)
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
A class for handling words, derived from Foam::string.
label size() const noexcept
The number of elements in the list.
constexpr scalar pi(M_PI)
A triFace with additional (region) index.
const Field< point_type > & points() const noexcept
Return reference to global points.
virtual const faceList & faces() const
Return raw faces.
bool typeHeaderOk(const bool checkType=true, const bool search=true, const bool verbose=true)
Read header (uses typeFilePath to find file) and check its info.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO...
const wordList surface
Standard surface field types (scalar, vector, tensor, etc)
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
static word timeName(const scalar t, const int precision=precision_)
Return time name of given scalar time formatted with the given precision.
const word & constant() const noexcept
Return constant name.
int debug
Static debugging option.
label nSolutionD() const
Return the number of valid solved-for dimensions in the mesh.
const boundBox & bounds() const
Return mesh bounding box.
label toGlobal(const label i) const
From local to global index.
Class containing processor-to-processor mapping information.
const word & name() const
Return reference to name.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A List with indirect addressing. Like IndirectList but does not store addressing. ...
static bool master(const label communicator=worldComm)
Am I the master rank.
label start() const
Return start label of this patch in the polyMesh face list.
const polyBoundaryMesh & patches
Automatically write from objectRegistry::writeObject()
const std::string patch
OpenFOAM patch number as a std::string.
vector span() const
The bounding box span (from minimum to maximum)
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Mesh consisting of general polyhedral cells.
A patch is a list of labels that address the faces in the global face list.
virtual bool write(const bool valid=true) const
Write using setting from DB.
Triangulated surface description with patch information.
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
Defines the attributes of an object for which implicit objectRegistry management is supported...
const Boundary & boundaryField() const noexcept
Return const-reference to the boundary field.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
forAllConstIters(mixture.phases(), phase)
CGAL::Segment_3< K > Segment
label localSize() const
My local size.
static constexpr const zero Zero
Global zero (0)