99 #include "CGAL/version.h" 100 #if defined(CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1050211000) 101 #define BOOST_BIND_GLOBAL_PLACEHOLDERS 103 #pragma clang diagnostic ignored "-Wbitwise-instead-of-logical" 106 #include <CGAL/Simple_cartesian.h> 107 #include <CGAL/AABB_tree.h> 108 #include <CGAL/AABB_traits.h> 109 #include <CGAL/AABB_triangle_primitive.h> 110 #include <CGAL/Surface_mesh.h> 112 typedef CGAL::Simple_cartesian<double>
K;
113 typedef K::Point_3
Point;
114 typedef K::Direction_3 Vector3C;
118 typedef std::vector<Triangle>::iterator Iterator;
119 typedef CGAL::AABB_triangle_primitive<K, Iterator> Primitive;
120 typedef CGAL::AABB_traits<K, Primitive> AABB_triangle_traits;
121 typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
122 typedef boost::optional
124 Tree::Intersection_and_primitive_id<Segment>::Type
125 > Segment_intersection;
129 using namespace Foam;
149 label localTriFaceI = 0;
151 for (
const label patchI : includePatches)
166 f.triangles(
points, nTri, triFaces);
168 forAll(triFaces, triFaceI)
170 const face&
f = triFaces[triFaceI];
179 finalAgglom[patchI][patchFaceI]
180 + coarsePatches[patchI].start()
214 for (
const label proci : globalPointIdx.allProcs())
216 const label offset = globalPointIdx.localStart(proci);
223 : globalSurfaceTris.slice(globalFaceIdx.range(proci))
236 std::move(globalSurfaceTris),
237 std::move(globalSurfacePoints)
248 for (
const label patchI : includePatches)
252 surface.patches()[newPatchI].index() = patchI;
254 surface.patches()[newPatchI].geometricType() =
patch.type();
278 const labelList visFaces = visibleFaceFaces[faceI];
279 forAll(visFaces, faceRemote)
281 label compactI = visFaces[faceRemote];
282 const point& remoteFc = compactCf[compactI];
288 str <<
"l " << vertI-1 <<
' ' << vertI <<
nl;
295 scalar calculateViewFactorFij2AI
304 scalar rMag =
mag(r);
308 scalar dAiMag =
mag(dAi);
309 scalar dAjMag =
mag(dAj);
313 scalar cosThetaJ =
mag(nj & r)/rMag;
314 scalar cosThetaI =
mag(ni & r)/rMag;
318 (cosThetaI*cosThetaJ*dAjMag*dAiMag)
329 void insertMatrixElements
332 const label fromProcI,
338 forAll(viewFactors, faceI)
341 const labelList& globalFaces = globalFaceFaces[faceI];
343 label globalI = globalNumbering.
toGlobal(fromProcI, faceI);
346 matrix[globalI][globalFaces[i]] = vf[i];
374 r = (
alpha*magSi)*di;
383 pi[i] = ci +
p[i]*(magSi/2)*di;
389 pj[i] = cj +
p[i]*(magSj/2)*dj;
400 r = (
alpha*magSi)*di;
413 dIntFij *= (magSi/2) * (magSj/2);
422 int main(
int argc,
char *argv[])
426 "Calculate view factors from face agglomeration array." 427 " The finalAgglom generated by faceAgglomerate utility." 448 const word viewFactorWall(
"viewFactorWall");
450 const bool writeViewFactors =
451 viewFactorDict.getOrDefault(
"writeViewFactorMatrix",
false);
453 const bool dumpRays =
454 viewFactorDict.getOrDefault(
"dumpRays",
false);
456 const label
debug = viewFactorDict.getOrDefault<label>(
"debug", 0);
458 const scalar GaussQuadTol =
459 viewFactorDict.getOrDefault<scalar>(
"GaussQuadTol", 0.01);
461 const scalar distTol =
462 viewFactorDict.getOrDefault<scalar>(
"distTol", 8);
467 const scalar intTol =
468 viewFactorDict.getOrDefault<scalar>(
"intTol", 1
e-2);
470 bool useAgglomeration(
true);
492 for (label patchi=0; patchi <
patches.
size(); patchi++)
496 useAgglomeration =
false;
504 Pout <<
"\nCreating single cell mesh..." <<
endl;
523 Pout <<
"\nCreated single cell mesh..." <<
endl;
530 label nCoarseFaces = 0;
531 label nFineFaces = 0;
535 for (
const label patchi : viewFactorsPatches)
537 nCoarseFaces += coarsePatches[patchi].
size();
542 Info<<
"\nTotal number of coarse faces: " 548 Pout <<
"\nView factor patches included in the calculation : " 549 << viewFactorsPatches <<
endl;
560 for (
const label
patchID : viewFactorsPatches)
567 if (agglom.
size() > 0)
569 label nAgglom =
max(agglom)+1;
572 coarseMesh.patchFaceMap()[
patchID];
575 coarseMesh.Cf().boundaryField()[
patchID];
577 coarseMesh.Sf().boundaryField()[
patchID];
581 point cf = coarseCf[faceI];
583 const label coarseFaceI = coarsePatchFace[faceI];
584 const labelList& fineFaces = coarseToFine[coarseFaceI];
585 const label agglomI =
597 upp.faceCentres().size()
598 + upp.localPoints().size()
604 upp.faceCentres().
size()
605 ) = upp.faceCentres();
610 upp.localPoints().size(),
611 upp.faceCentres().size()
612 ) = upp.localPoints();
616 forAll(availablePoints, iPoint)
618 point cfFine = availablePoints[iPoint];
619 if (
mag(cfFine-cfo) < dist)
621 dist =
mag(cfFine-cfo);
626 point sf = coarseSf[faceI];
628 localCoarseSf.append(sf);
629 localAgg.append(agglomI);
684 nVisibleFaceFaces[rayStartFace[i]]++;
689 label nViewFactors = 0;
690 forAll(nVisibleFaceFaces, faceI)
692 visibleFaceFaces[faceI].
setSize(nVisibleFaceFaces[faceI]);
693 nViewFactors += nVisibleFaceFaces[faceI];
710 nVisibleFaceFaces = 0;
713 label faceI = rayStartFace[i];
714 label compactI = rayEndFace[i];
715 visibleFaceFaces[faceI][nVisibleFaceFaces[faceI]++] = compactI;
740 forAll(viewFactorsPatches, i)
742 label
patchID = viewFactorsPatches[i];
745 if (agglom.
size() > 0)
747 label nAgglom =
max(agglom)+1;
750 coarseMesh.patchFaceMap()[
patchID];
754 forAll(coarseToFine, coarseI)
756 compactPatchId.append(
patchID);
760 label startFace =
pp.start();
765 faces[coarseI + startFace]
768 const label coarseFaceI = coarsePatchFace[coarseI];
769 const labelList& fineFaces = coarseToFine[coarseFaceI];
774 compactPoints.append(locPoints);
779 coarseToFine[coarseFaceI]
784 coarseToFine[coarseFaceI]
798 map.distribute(compactCoarseSf);
799 map.distribute(compactCoarseCf);
800 map.distribute(compactFineCf);
801 map.distribute(compactFineSf);
802 map.distribute(compactPoints);
803 map.distribute(compactPatchId);
834 const label totalPatches =
844 Info<<
"\nCalculating view factors..." <<
endl;
856 GaussPoints[2][0] = 0;
862 GaussPoints[3][1] = -GaussPoints[3][0];
864 GaussPoints[3][3] = -GaussPoints[3][2];
867 GaussPoints[4][0] = 0;
869 GaussPoints[4][2] = -GaussPoints[4][1];
871 GaussPoints[4][4] = -GaussPoints[4][3];
879 GaussWeights[1][0] = 1;
880 GaussWeights[1][1] = 1;
883 GaussWeights[2][0] = 8.0/9.0;
884 GaussWeights[2][1] = 5.0/9.0;
885 GaussWeights[2][2] = 5.0/9.0;
888 GaussWeights[3][0] = (18.0 +
std::sqrt(30))/36.0;
889 GaussWeights[3][1] = (18.0 +
std::sqrt(30))/36.0;
890 GaussWeights[3][2] = (18.0 -
std::sqrt(30))/36.0;
891 GaussWeights[3][3] = (18.0 -
std::sqrt(30))/36.0;
894 GaussWeights[4][0] = 128.0/225.0;
895 GaussWeights[4][1] = (322.0 + 13.0*
std::sqrt(70))/900.0;
896 GaussWeights[4][2] = (322.0 + 13.0*
std::sqrt(70))/900.0;
897 GaussWeights[4][3] = (322.0 - 13.0*
std::sqrt(70))/900.0;
898 GaussWeights[4][4] = (322.0 - 13.0*
std::sqrt(70))/900.0;
900 const label maxQuadOrder = 5;
904 forAll(localCoarseSf, coarseFaceI)
906 const List<point>& localFineSf = compactFineSf[coarseFaceI];
908 const List<point>& localFineCf = compactFineCf[coarseFaceI];
909 const label fromPatchId = compactPatchId[coarseFaceI];
911 const List<point>& lPoints = compactPoints[coarseFaceI];
913 patchArea[fromPatchId] +=
mag(Ai);
915 const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
917 forAll(visCoarseFaces, visCoarseFaceI)
920 F2LI[coarseFaceI].setSize(visCoarseFaces.
size());
921 label compactJ = visCoarseFaces[visCoarseFaceI];
922 const List<point>& remoteFineSj = compactFineSf[compactJ];
923 const List<point>& remoteFineCj = compactFineCf[compactJ];
925 const List<point>& rPointsCj = compactPoints[compactJ];
927 const label toPatchId = compactPatchId[compactJ];
938 const vector& dCi = localFineCf[i];
947 const vector& dCj = remoteFineCj[j];
949 const scalar dist =
mag(dCi - dCj)/((dAi + dAj)/2);
965 const vector& dAi = localFineSf[i];
966 const vector& dCi = localFineCf[i];
970 const vector& dAj = remoteFineSj[j];
971 const vector& dCj = remoteFineCj[j];
973 scalar dIntFij = calculateViewFactorFij2AI
984 F2LI[coarseFaceI][visCoarseFaceI] = F2AIij/
mag(Ai);
989 label nLocal = lPoints.
size();
990 label nRemote = rPointsCj.
size();
993 scalar oldEToeInt = 0;
994 for (label gi=0; gi < maxQuadOrder; gi++)
997 for(label i=0; i<nLocal; i++)
1007 si = lPoints[i] - lPoints[nLocal-1];
1008 ci = (lPoints[i] + lPoints[nLocal-1])/2;
1012 si = lPoints[i] - lPoints[i-1];
1013 ci = (lPoints[i] + lPoints[i-1])/2;
1016 for(label j=0; j<nRemote; j++)
1020 sj = rPointsCj[j]-rPointsCj[nRemote-1];
1021 cj = (rPointsCj[j]+rPointsCj[nRemote-1])/2;
1025 sj = rPointsCj[j] - rPointsCj[j-1];
1026 cj = (rPointsCj[j] + rPointsCj[j-1])/2;
1030 scalar magSi =
mag(si);
1031 scalar magSj =
mag(sj);
1032 scalar cosij = (si & sj)/(magSi * magSj);
1037 label quadOrder = gi;
1067 mag(F2LIij) > ROOTVSMALL
1068 ? (F2LIij-oldEToeInt)/F2LIij
1074 (
mag(err) < GaussQuadTol && gi > 0)
1075 || gi == maxQuadOrder-1
1078 F2LI[coarseFaceI][visCoarseFaceI] =
1084 oldEToeInt = F2LIij;
1089 sumViewFactorPatch[fromPatchId][toPatchId] +=
1090 F2LI[coarseFaceI][visCoarseFaceI]*
mag(Ai);
1107 scalar wideBy2 = (box.
span() & emptyDir)*2.0;
1109 forAll(localCoarseSf, coarseFaceI)
1111 const vector& Ai = localCoarseSf[coarseFaceI];
1112 const vector& Ci = localCoarseCf[coarseFaceI];
1114 vector R1i = Ci + (
mag(Ai)/wideBy2)*(Ain ^ emptyDir);
1115 vector R2i = Ci - (
mag(Ai)/wideBy2)*(Ain ^ emptyDir) ;
1117 const label fromPatchId = compactPatchId[coarseFaceI];
1118 patchArea[fromPatchId] +=
mag(Ai);
1120 const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
1121 forAll(visCoarseFaces, visCoarseFaceI)
1123 F2LI[coarseFaceI].setSize(visCoarseFaces.
size());
1124 label compactJ = visCoarseFaces[visCoarseFaceI];
1125 const vector& Aj = compactCoarseSf[compactJ];
1126 const vector& Cj = compactCoarseCf[compactJ];
1128 const label toPatchId = compactPatchId[compactJ];
1131 vector R1j = Cj + (
mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
1132 vector R2j = Cj - (
mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
1134 scalar d1 =
mag(R1i - R2j);
1135 scalar d2 =
mag(R2i - R1j);
1136 scalar s1 =
mag(R1i - R1j);
1137 scalar s2 =
mag(R2i - R2j);
1139 scalar Fij =
mag((d1 + d2) - (s1 + s2))/(4.0*
mag(Ai)/wideBy2);
1141 F2LI[coarseFaceI][visCoarseFaceI] = Fij;
1142 sumViewFactorPatch[fromPatchId][toPatchId] += Fij*
mag(Ai);
1149 <<
" View factors are not available in 1D " 1161 forAll(viewFactorsPatches, i)
1163 label patchI = viewFactorsPatches[i];
1164 for (label j=i; j<viewFactorsPatches.size(); j++)
1166 label patchJ = viewFactorsPatches[j];
1168 Info <<
"F" << patchI << patchJ <<
": " 1169 << sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
1175 if (writeViewFactors)
1179 Info <<
"Writing view factor matrix..." <<
endl;
1199 forAll(viewFactorsPatches, i)
1201 label
patchID = viewFactorsPatches[i];
1203 if (agglom.
size() > 0)
1205 label nAgglom =
max(agglom)+1;
1208 coarseMesh.patchFaceMap()[
patchID];
1210 forAll(coarseToFine, coarseI)
1212 const scalar FiSum =
sum(F2LI[compactI]);
1214 const label coarseFaceID = coarsePatchFace[coarseI];
1215 const labelList& fineFaces = coarseToFine[coarseFaceID];
1216 forAll(fineFaces, fineId)
1218 const label faceID = fineFaces[fineId];
1219 vfbf[
patchID][faceID] = FiSum;
1225 viewFactorField.write();
1232 labelList compactToGlobal(map.constructSize());
1235 for (label i = 0; i < globalNumbering.
localSize(); i++)
1237 compactToGlobal[i] = globalNumbering.
toGlobal(i);
1241 forAll(compactMap, procI)
1243 const Map<label>& localToCompactMap = compactMap[procI];
1247 compactToGlobal[*iter] = globalNumbering.
toGlobal 1260 forAll(globalFaceFaces, faceI)
1265 visibleFaceFaces[faceI]
1280 std::move(globalFaceFaces)
1282 IOglobalFaceFaces.
write();
Different types of constants.
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.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
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 = rootPath/caseName. Same as TimePaths::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)
Rank of this process in the communicator (starting from masterNo()). Can be negative if the process i...
const Time & time() const
Return the top-level database.
virtual const fileName & name() const override
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
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 communicator ranks. Does nothing in non-paral...
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). It 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), works like std::iota() but returning a...
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
A class for handling words, derived from Foam::string.
label size() const noexcept
The number of entries in the list.
Reading is optional [identical to LAZY_READ].
virtual bool write(const bool writeOnProc=true) const
Write using setting from DB.
constexpr scalar pi(M_PI)
label localSize(const label proci) const
Size of proci data.
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 (respects is_globalIOobject trait) 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 a time name for the given scalar time value 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.
Class containing processor-to-processor mapping information.
label toGlobal(const label proci, const label i) const
From local to global on proci.
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)
True if process corresponds to the master rank in the communicator.
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)
const boundBox & bounds() const noexcept
Return mesh bounding box.
Mesh consisting of general polyhedral cells.
A patch is a list of labels that address the faces in the global face list.
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.
Do not request registration (bool: false)
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
forAllConstIters(mixture.phases(), phase)
CGAL::Segment_3< K > Segment
static constexpr const zero Zero
Global zero (0)