diff --git a/applications/utilities/preProcessing/viewFactorsGen/Make/options b/applications/utilities/preProcessing/viewFactorsGen/Make/options
index 1a9826a229ce5cbe4b55e2ccd1bda80fcbeec3c2..298501f44781b272f32ae13c91f675c9498ae9e6 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/Make/options
+++ b/applications/utilities/preProcessing/viewFactorsGen/Make/options
@@ -1,10 +1,16 @@
+include $(GENERAL_RULES)/cgal
+
EXE_INC = \
+ -Wno-old-style-cast \
+ ${CGAL_INC} \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
- -I$(LIB_SRC)/parallel/distributed/lnInclude \
+ -I$(LIB_SRC)/parallel/distributed/lnInclude
+
EXE_LIBS = \
+ ${CGAL_LIBS} \
-lfiniteVolume \
-lsurfMesh \
-lmeshTools \
diff --git a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
index df38326ded081fc52a391c42f5cf01212e5e3568..9b89b848c989f13953207c5ff1a611b9ce42eba8 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
@@ -25,7 +25,7 @@ dict.add("mergeDistance", SMALL);
labelList triSurfaceToAgglom(5*nFineFaces);
-const triSurface localSurface = triangulate
+triSurface localSurface = triangulate
(
patches,
includePatches,
@@ -36,6 +36,8 @@ const triSurface localSurface = triangulate
);
+// CGAL surface
+
distributedTriSurfaceMesh surfacesMesh
(
IOobject
@@ -45,13 +47,38 @@ distributedTriSurfaceMesh surfacesMesh
"triSurface", // instance
runTime, // registry
IOobject::NO_READ,
- IOobject::NO_WRITE
+ IOobject::AUTO_WRITE
),
localSurface,
dict
);
-
triSurfaceToAgglom.resize(surfacesMesh.size());
surfacesMesh.setField(triSurfaceToAgglom);
+
+const pointField& pts = surfacesMesh.localPoints();
+
+std::list<Triangle> triangles;
+
+for (const auto& triLocal : surfacesMesh.localFaces())
+{
+ point p1l = pts[triLocal[0]];
+ point p2l = pts[triLocal[1]];
+ point p3l = pts[triLocal[2]];
+
+ Point p1(p1l[0], p1l[1], p1l[2]);
+ Point p2(p2l[0], p2l[1], p2l[2]);
+ Point p3(p3l[0], p3l[1], p3l[2]);
+
+ Triangle tri(p1, p2, p3);
+
+ if (tri.is_degenerate())
+ {
+ std::cout << tri << std::endl;
+ }
+ triangles.push_back(tri);
+}
+
+// constructs AABB tree
+Tree tree(triangles.begin(), triangles.end());
diff --git a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
index e35b2d4268c6ae08ac4be3e7e1ca68fe1af78dac..c738784b79939f6ee2554eabf9f7772c62df0819 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
@@ -1,26 +1,23 @@
-// All rays expressed as start face (local) index end end face (global)
-// Pre-size by assuming a certain percentage is visible.
-
// Maximum length for dynamicList
const label maxDynListLength
(
- viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000)
+ viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000000)
);
for (const int proci : Pstream::allProcs())
{
- // Shoot rays from me to proci. Note that even if processor has
- // 0 faces we still need to call findLine to keep calls synced.
+ std::vector<Point> start;
+ start.reserve(coarseMesh.nFaces());
+
+ std::vector<Point> end(start.size());
+ end.reserve(start.size());
- DynamicField<point> start(coarseMesh.nFaces());
- DynamicField<point> end(start.size());
DynamicList<label> startIndex(start.size());
DynamicList<label> endIndex(start.size());
DynamicList<label> startAgg(start.size());
DynamicList<label> endAgg(start.size());
-
const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
@@ -50,15 +47,30 @@ for (const int proci : Pstream::allProcs())
const vector d(remFc - fc);
- if (((d & fA) < 0.) && ((d & remA) > 0))
+ const vector nd = d/mag(d);
+ const vector nfA = fA/mag(fA);
+ const vector nremA = remA/mag(remA);
+
+ if (((nd & nfA) < 0) && ((nd & nremA) > 0))
{
- start.append(fc + 0.001*d);
+ vector direction(d[0], d[1], d[2]);
+ vector s(fc[0], fc[1], fc[2]);
+ vector rayEnd(s + (1-intTol)*direction);
+ end.push_back(Point(rayEnd[0], rayEnd[1], rayEnd[2]));
+
+ s += vector(intTol*d[0], intTol*d[1], intTol*d[2]);
+
+ start.push_back(Point(s[0], s[1], s[2]));
startIndex.append(i);
- startAgg.append(globalNumbering.toGlobal(proci, fAgg));
- end.append(fc + 0.999*d);
+ if (useAgglomeration)
+ {
+ startAgg.append(globalNumbering.toGlobal(proci, fAgg));
+ endAgg.append(globalNumbering.toGlobal(proci, remAgg));
+ }
+
label globalI = globalNumbering.toGlobal(proci, j);
endIndex.append(globalI);
- endAgg.append(globalNumbering.toGlobal(proci, remAgg));
+
if (startIndex.size() > maxDynListLength)
{
FatalErrorInFunction
@@ -77,98 +89,33 @@ for (const int proci : Pstream::allProcs())
}
}
- } while (returnReduce(i < myFc.size(), orOp<bool>()));
-
- List<pointIndexHit> hitInfo(startIndex.size());
- surfacesMesh.findLine(start, end, hitInfo);
+ } while (i < myFc.size());
- // Return hit global agglo index
- labelList aggHitIndex;
- surfacesMesh.getField(hitInfo, aggHitIndex);
+ label totalRays(startIndex.size());
+ reduce(totalRays, sumOp<scalar>());
- DynamicList<label> dRayIs;
+ if (debug)
+ {
+ Pout<< "Number of rays :" << totalRays << endl;
+ }
- // Collect the rays which have no obstacle in between rayStartFace
- // and rayEndFace. If the ray hit itself, it gets stored in dRayIs
- forAll(hitInfo, rayI)
+ for (unsigned long rayI = 0; rayI < start.size(); ++rayI)
{
- if (!hitInfo[rayI].hit())
+ Segment ray(start[rayI], end[rayI]);
+ Segment_intersection intersects = tree.any_intersection(ray);
+
+ if (!intersects)
{
rayStartFace.append(startIndex[rayI]);
rayEndFace.append(endIndex[rayI]);
}
- else if (aggHitIndex[rayI] == startAgg[rayI])
- {
- dRayIs.append(rayI);
- }
}
-
start.clear();
-
- // Continue rays which hit themself. If they hit the target
- // agglomeration are added to rayStartFace and rayEndFace
-
- bool firstLoop = true;
- DynamicField<point> startHitItself;
- DynamicField<point> endHitItself;
- label iter = 0;
-
- do
+ if (debug)
{
- labelField rayIs;
- rayIs.transfer(dRayIs);
- dRayIs.clear();
- forAll(rayIs, rayI)
- {
- const label rayID = rayIs[rayI];
- label hitIndex = -1;
-
- if (firstLoop)
- {
- hitIndex = rayIs[rayI];
- }
- else
- {
- hitIndex = rayI;
- }
-
- if (hitInfo[hitIndex].hit())
- {
- if (aggHitIndex[hitIndex] == startAgg[rayID])
- {
- const vector& endP = end[rayID];
- const vector& startP = hitInfo[hitIndex].hitPoint();
- const vector& d = endP - startP;
-
- startHitItself.append(startP + 0.01*d);
- endHitItself.append(startP + 1.01*d);
-
- dRayIs.append(rayID);
- }
- else if (aggHitIndex[hitIndex] == endAgg[rayID])
- {
- rayStartFace.append(startIndex[rayID]);
- rayEndFace.append(endIndex[rayID]);
- }
-
- }
- }
-
- hitInfo.clear();
- hitInfo.resize(dRayIs.size());
-
- surfacesMesh.findLine(startHitItself, endHitItself, hitInfo);
-
- surfacesMesh.getField(hitInfo, aggHitIndex);
-
-
- endHitItself.clear();
- startHitItself.clear();
- firstLoop = false;
- iter ++;
-
- } while (returnReduce(hitInfo.size(), orOp<bool>()) > 0 && iter < 10);
+ Pout << "hits : " << rayStartFace.size()<< endl;
+ }
startIndex.clear();
end.clear();
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index fafdfb811b9d3cf9e757efefdf897c8bacf9fbee..368fc8df1c78cf39383ad8ba8abae2ce7f30b628 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -4,7 +4,7 @@
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
--------------------------------------------------------------------------------
+// -------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2016-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
@@ -25,20 +25,40 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
- viewFactorsGen
+ viewFactorsGenExt
Group
grpPreProcessingUtilities
Description
- View factors are calculated based on a face agglomeration array
- (finalAgglom generated by faceAgglomerate utility).
+ This view factors generation application uses a combined approach of
+ double area integral (2AI) and double linear integral (2LI). 2AI is used
+ when the two surfaces are 'far' apart and 2LI whenre they are 'close'.
+ 2LI is integrated along edges using Gaussian quadrature.
+ The distance between faces is calculating a ratio between averaged areas
+ and the distance between face centres.
+
+ The input from viewFactorsDict are:
+
+ tolGaussQuad 0.1; // GaussQuad error
+ distTol 8; // R/Average(rm)
+ alpha 0.22; // Use for common edges for 2LI
+
+
+ For debugging purposes, the following entries can be set in viewFactorsDict:
- Each view factor between the agglomerated faces i and j (Fij) is calculated
- using a double integral of the sub-areas composing the agglomeration.
+ writeViewFactorMatrix true;
+ writeFacesAgglomeration false;
+ dumpRays false;
- The patches involved in the view factor calculation are taken from the
- boundary file and should be part on the group viewFactorWall. ie.:
+
+ writeViewFactorMatrix writes the sum of the VF on each face.
+ writeFacesAgglomeration writes the agglomeration
+ dumpRays dumps rays
+
+
+ In order to specify the participants patches in the VF calculation the
+ keywaord viewFactorWall should be added to the boundary file.
floor
{
@@ -56,21 +76,55 @@ Description
#include "surfaceFields.H"
#include "distributedTriSurfaceMesh.H"
#include "meshTools.H"
+#include "constants.H"
#include "indirectPrimitivePatch.H"
#include "DynamicField.H"
-#include "unitConversion.H"
+//#include "unitConversion.H"
#include "scalarMatrices.H"
#include "labelListIOList.H"
#include "scalarListIOList.H"
#include "singleCellFvMesh.H"
-
#include "IOmapDistribute.H"
+#ifndef NO_CGAL
+
+// Silence boost bind deprecation warnings (before CGAL-5.2.1)
+#include "CGAL/version.h"
+#if defined(CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1050211000)
+#define BOOST_BIND_GLOBAL_PLACEHOLDERS
+#endif
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/Surface_mesh.h>
+
+typedef CGAL::Simple_cartesian<double> K;
+typedef K::Point_3 Point;
+typedef K::Direction_3 Vector3C;
+typedef K::Triangle_3 Triangle;
+typedef K::Segment_3 Segment;
+
+typedef std::list<Triangle>::iterator Iterator;
+typedef CGAL::AABB_triangle_primitive<K, Iterator> Primitive;
+typedef CGAL::AABB_traits<K, Primitive> AABB_triangle_traits;
+typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+typedef boost::optional
+<
+ Tree::Intersection_and_primitive_id<Segment>::Type
+> Segment_intersection;
+
+#endif // NO_CGAL
+
using namespace Foam;
+using namespace Foam::constant;
+using namespace Foam::constant::mathematical;
+//using namespace pbrt;
triSurface triangulate
(
@@ -127,20 +181,74 @@ triSurface triangulate
newPatchI++;
}
- //striSurfaceToAgglom.resize(localTriFaceI-1);
+ //triSurfaceToAgglom.resize(localTriFaceI-1);
triangles.shrink();
- // Create globally numbered tri surface
triSurface rawSurface(triangles, mesh.points());
- // Create locally numbered tri surface
triSurface surface
(
rawSurface.localFaces(),
rawSurface.localPoints()
);
+ // Combine the triSurfaces across all processors
+ if (Pstream::parRun())
+ {
+ List<List<labelledTri>> surfaceProcTris(Pstream::nProcs());
+ List<pointField> surfaceProcPoints(Pstream::nProcs());
+
+ surfaceProcTris[Pstream::myProcNo()] = surface;
+ surfaceProcPoints[Pstream::myProcNo()] = surface.points();
+
+ Pstream::gatherList(surfaceProcTris);
+ Pstream::scatterList(surfaceProcTris);
+ Pstream::gatherList(surfaceProcPoints);
+ Pstream::scatterList(surfaceProcPoints);
+
+ label nTris = 0;
+ forAll(surfaceProcTris, i)
+ {
+ nTris += surfaceProcTris[i].size();
+ }
+
+ List<labelledTri> globalSurfaceTris(nTris);
+ label trii = 0;
+ label offset = 0;
+ forAll(surfaceProcTris, i)
+ {
+ forAll(surfaceProcTris[i], j)
+ {
+ globalSurfaceTris[trii] = surfaceProcTris[i][j];
+ globalSurfaceTris[trii][0] += offset;
+ globalSurfaceTris[trii][1] += offset;
+ globalSurfaceTris[trii][2] += offset;
+ trii++;
+ }
+ offset += surfaceProcPoints[i].size();
+ }
+
+ label nPoints = 0;
+ forAll(surfaceProcPoints, i)
+ {
+ nPoints += surfaceProcPoints[i].size();
+ }
+
+ pointField globalSurfacePoints(nPoints);
+
+ label pointi = 0;
+ forAll(surfaceProcPoints, i)
+ {
+ forAll(surfaceProcPoints[i], j)
+ {
+ globalSurfacePoints[pointi++] = surfaceProcPoints[i][j];
+ }
+ }
+
+ surface = triSurface(globalSurfaceTris, globalSurfacePoints);
+ }
+
// Add patch names to surface
surface.patches().setSize(newPatchI);
@@ -193,7 +301,7 @@ void writeRays
}
-scalar calculateViewFactorFij
+scalar calculateViewFactorFij2AI
(
const vector& i,
const vector& j,
@@ -249,7 +357,72 @@ void insertMatrixElements
}
}
+scalar GaussQuad
+(
+ const scalarList& w,
+ const scalarList& p,
+ const scalar& magSi,
+ const scalar& magSj,
+ const vector& di,
+ const vector& dj,
+ const vector& ci,
+ const vector& cj,
+ const scalar cosij,
+ const scalar alpha,
+ label gi
+)
+{
+ scalar dIntFij = 0;
+ if (gi == 0)
+ {
+ vector r(ci - cj);
+ if (mag(r) < SMALL)
+ {
+ r = (alpha*magSi)*di;
+ }
+ dIntFij = max(cosij*Foam::log(r&r)*magSi*magSj, 0);
+ }
+ else
+ {
+ List<vector> pi(w.size());
+ forAll (pi, i)
+ {
+ pi[i] = ci + p[i]*(magSi/2)*di;
+ }
+
+ List<vector> pj(w.size());
+ forAll (pj, i)
+ {
+ pj[i] = cj + p[i]*(magSj/2)*dj;
+ }
+
+ forAll (w, i)
+ {
+ forAll (w, j)
+ {
+ vector r(pi[i] - pj[j]);
+
+ if (mag(r) < SMALL)
+ {
+ r = (alpha*magSi)*di;
+ dIntFij +=
+ cosij*w[i]*w[j]*Foam::log(r&r);
+ }
+ else
+ {
+ dIntFij +=
+ cosij*w[i]*w[j]*Foam::log(r&r);
+ }
+
+ }
+ }
+
+ dIntFij *= (magSi/2) * (magSj/2);
+
+ }
+ return dIntFij;
+}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
@@ -273,7 +446,7 @@ int main(int argc, char *argv[])
"viewFactorsDict",
runTime.constant(),
mesh,
- IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
@@ -288,6 +461,23 @@ int main(int argc, char *argv[])
const label debug = viewFactorDict.getOrDefault<label>("debug", 0);
+ const scalar tolGaussQuad =
+ viewFactorDict.getOrDefault<scalar>("tolGaussQuad", 0.01);
+
+ const scalar distTol =
+ viewFactorDict.getOrDefault<scalar>("distTol", 8);
+
+ const scalar alpha =
+ viewFactorDict.getOrDefault<scalar>("alpha", 0.21);
+
+ const scalar intTol =
+ viewFactorDict.getOrDefault<scalar>("intTol", 1e-2);
+
+ bool useAgglomeration(true);
+
+ const polyBoundaryMesh& patches = mesh.boundaryMesh();
+ const labelList viewFactorsPatches(patches.indices(viewFactorWall));
+
// Read agglomeration map
labelListIOList finalAgglom
(
@@ -296,12 +486,22 @@ int main(int argc, char *argv[])
"finalAgglom",
mesh.facesInstance(),
mesh,
- IOobject::MUST_READ,
+ IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE,
false
)
);
+ if (!finalAgglom.typeHeaderOk<labelListIOList>())
+ {
+ finalAgglom.setSize(patches.size());
+ for (label patchi=0; patchi < patches.size(); patchi++)
+ {
+ finalAgglom[patchi] = identity(patches[patchi].size());
+ }
+ useAgglomeration = false;
+ }
+
// Create the coarse mesh using agglomeration
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -336,10 +536,8 @@ int main(int argc, char *argv[])
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.indices(viewFactorWall));
for (const label patchi : viewFactorsPatches)
{
nCoarseFaces += coarsePatches[patchi].size();
@@ -473,10 +671,8 @@ int main(int argc, char *argv[])
DynamicList<label> rayEndFace(rayStartFace.size());
-
// Return rayStartFace in local index and rayEndFace in global index
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
#include "shootRays.H"
// Calculate number of visible faces from local index
@@ -518,9 +714,9 @@ int main(int argc, char *argv[])
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
+ // (2AA) need coarse (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
+ // (2LI) need edges (li)
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
pointField compactCoarseCf(map.constructSize(), Zero);
@@ -528,12 +724,16 @@ int main(int argc, char *argv[])
List<List<point>> compactFineSf(map.constructSize());
List<List<point>> compactFineCf(map.constructSize());
+ DynamicList<List<point>> compactPoints(map.constructSize());
+
DynamicList<label> compactPatchId(map.constructSize());
// Insert my coarse local values
SubList<point>(compactCoarseSf, nCoarseFaces) = localCoarseSf;
SubList<point>(compactCoarseCf, nCoarseFaces) = localCoarseCf;
+ const faceList& faces = mesh.faces();
+
// Insert my fine local values
label compactI = 0;
forAll(viewFactorsPatches, i)
@@ -548,11 +748,21 @@ int main(int argc, char *argv[])
const labelList& coarsePatchFace =
coarseMesh.patchFaceMap()[patchID];
+ const polyPatch& pp = patches[patchID];
+
forAll(coarseToFine, coarseI)
{
compactPatchId.append(patchID);
List<point>& fineCf = compactFineCf[compactI];
- List<point>& fineSf = compactFineSf[compactI++];
+ List<point>& fineSf = compactFineSf[compactI];
+
+ label startFace = pp.start();
+
+ const vectorField locPoints
+ (
+ mesh.points(),
+ faces[coarseI + startFace]
+ );
const label coarseFaceI = coarsePatchFace[coarseI];
const labelList& fineFaces = coarseToFine[coarseFaceI];
@@ -560,6 +770,8 @@ int main(int argc, char *argv[])
fineCf.setSize(fineFaces.size());
fineSf.setSize(fineFaces.size());
+ compactPoints.append(locPoints);
+
fineCf = UIndirectList<point>
(
mesh.Cf().boundaryField()[patchID],
@@ -570,19 +782,25 @@ int main(int argc, char *argv[])
mesh.Sf().boundaryField()[patchID],
coarseToFine[coarseFaceI]
);
+
+ compactI++;
}
}
}
+ if (Pstream::master() && debug)
+ {
+ Info<< "map distribute..." << endl;
+ }
+
// Do all swapping
map.distribute(compactCoarseSf);
map.distribute(compactCoarseCf);
map.distribute(compactFineCf);
map.distribute(compactFineSf);
-
+ map.distribute(compactPoints);
map.distribute(compactPatchId);
-
// Plot all rays between visible faces.
if (dumpRays)
{
@@ -598,8 +816,7 @@ int main(int argc, char *argv[])
// Fill local view factor matrix
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- scalarListIOList F
+ scalarListIOList F2LI
(
IOobject
(
@@ -630,6 +847,61 @@ int main(int argc, char *argv[])
Info<< "\nCalculating view factors..." << endl;
}
+ FixedList<scalarList, 5> GaussPoints;
+ GaussPoints[0].setSize(1);
+ GaussPoints[0] = 0;
+
+ GaussPoints[1].setSize(2);
+ GaussPoints[1][0] = 1/std::sqrt(3);
+ GaussPoints[1][1] = -1/std::sqrt(3);
+
+ GaussPoints[2].setSize(3);
+ GaussPoints[2][0] = 0;
+ GaussPoints[2][1] = 0.774597;
+ GaussPoints[2][2] = -0.774597;
+
+ GaussPoints[3].setSize(4);
+ GaussPoints[3][0] = 0.339981;
+ GaussPoints[3][1] = -0.339981;
+ GaussPoints[3][2] = 0.861136;
+ GaussPoints[3][3] = -0.861136;
+
+ GaussPoints[4].setSize(5);
+ GaussPoints[4][0] = 0;
+ GaussPoints[4][1] = 0.538469;
+ GaussPoints[4][2] = -0.538469;
+ GaussPoints[4][3] = 0.90618;
+ GaussPoints[4][4] = -0.90618;
+
+
+ List<scalarList> GaussWeights(5);
+ GaussWeights[0].setSize(1);
+ GaussWeights[0] = 2;
+
+ GaussWeights[1].setSize(2);
+ GaussWeights[1][0] = 1;
+ GaussWeights[1][1] = 1;
+
+ GaussWeights[2].setSize(3);
+ GaussWeights[2][0] = 0.888889;
+ GaussWeights[2][1] = 0.555556;
+ GaussWeights[2][2] = 0.555556;
+
+ GaussWeights[3].setSize(4);
+ GaussWeights[3][0] = 0.652145;
+ GaussWeights[3][1] = 0.652145;
+ GaussWeights[3][2] = 0.347855;
+ GaussWeights[3][3] = 0.347855;
+
+ GaussWeights[4].setSize(5);
+ GaussWeights[4][0] = 0.568889;
+ GaussWeights[4][1] = 0.478629;
+ GaussWeights[4][2] = 0.478629;
+ GaussWeights[4][3] = 0.236927;
+ GaussWeights[4][4] = 0.236927;
+
+ label maxQuadOrder = 5;
+
if (mesh.nSolutionD() == 3)
{
forAll(localCoarseSf, coarseFaceI)
@@ -638,119 +910,207 @@ int main(int argc, char *argv[])
const vector Ai = sum(localFineSf);
const List<point>& localFineCf = compactFineCf[coarseFaceI];
const label fromPatchId = compactPatchId[coarseFaceI];
+
+ const List<point>& lPoints = compactPoints[coarseFaceI];
+
patchArea[fromPatchId] += mag(Ai);
const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
forAll(visCoarseFaces, visCoarseFaceI)
{
- F[coarseFaceI].setSize(visCoarseFaces.size());
+ //F2AI[coarseFaceI].setSize(visCoarseFaces.size());
+ F2LI[coarseFaceI].setSize(visCoarseFaces.size());
label compactJ = visCoarseFaces[visCoarseFaceI];
const List<point>& remoteFineSj = compactFineSf[compactJ];
const List<point>& remoteFineCj = compactFineCf[compactJ];
+ const List<point>& rPointsCj = compactPoints[compactJ];
+
const label toPatchId = compactPatchId[compactJ];
- scalar Fij = 0;
+ bool far(false);
+ // Relative distance
forAll(localFineSf, i)
{
- const vector& dAi = localFineSf[i];
+ const scalar dAi =
+ Foam::sqrt
+ (
+ mag(localFineSf[i])/constant::mathematical::pi
+ );
const vector& dCi = localFineCf[i];
forAll(remoteFineSj, j)
{
- const vector& dAj = remoteFineSj[j];
+ const scalar dAj =
+ Foam::sqrt
+ (
+ mag(remoteFineSj[j])/constant::mathematical::pi
+ );
const vector& dCj = remoteFineCj[j];
- scalar dIntFij = calculateViewFactorFij
- (
- dCi,
- dCj,
- dAi,
- dAj
- );
+ const scalar dist = mag(dCi - dCj)/((dAi + dAj)/2);
- Fij += dIntFij;
+ if (dist > distTol)
+ {
+ far = true;
+ }
}
}
- 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 = 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);
+ if (far)
+ {
+ // 2AI method
+ scalar F2AIij = 0;
- scalar d1 = mag(R1i - R2j);
- scalar d2 = mag(R2i - R1j);
- scalar s1 = mag(R1i - R1j);
- scalar s2 = mag(R2i - R2j);
+ 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 = calculateViewFactorFij2AI
+ (
+ dCi,
+ dCj,
+ dAi,
+ dAj
+ );
+
+ F2AIij += dIntFij;
+ }
+ }
+ F2LI[coarseFaceI][visCoarseFaceI] = F2AIij/mag(Ai);
+ }
+ else
+ {
+ // 2LI method
+ label nLocal = lPoints.size();
+ label nRemote = rPointsCj.size();
- scalar Fij = mag((d1 + d2) - (s1 + s2))/(4.0*mag(Ai)/wideBy2);
+ // Using sub-divisions (quadrature)
+ scalar oldEToeInt = 0;
+ for (label gi=0; gi < maxQuadOrder; gi++)
+ {
+ scalar F2LIij = 0;
+ for(label i=0; i<nLocal; i++)
+ {
+ vector si;
+ vector ci;
+
+ vector sj;
+ vector cj;
+
+ if (i == 0)
+ {
+ si = lPoints[i] - lPoints[nLocal-1];
+ ci = (lPoints[i] + lPoints[nLocal-1])/2;
+ }
+ else
+ {
+ si = lPoints[i] - lPoints[i-1];
+ ci = (lPoints[i] + lPoints[i-1])/2;
+ }
+
+ for(label j=0; j<nRemote; j++)
+ {
+ if (j == 0)
+ {
+ sj = rPointsCj[j]-rPointsCj[nRemote-1];
+ cj = (rPointsCj[j]+rPointsCj[nRemote-1])/2;
+ }
+ else
+ {
+ sj = rPointsCj[j] - rPointsCj[j-1];
+ cj = (rPointsCj[j] + rPointsCj[j-1])/2;
+ }
+
+
+ scalar magSi = mag(si);
+ scalar magSj = mag(sj);
+ scalar cosij = (si & sj)/(magSi * magSj);
+
+ vector di = si/magSi;
+ vector dj = sj/magSj;
+
+ label quadOrder = gi;
+ const vector r(ci - cj);
+ // Common edges use n = 0
+ if (mag(r) < SMALL)
+ {
+ quadOrder = 0;
+ }
+
+ scalar dIntFij =
+ GaussQuad
+ (
+ GaussWeights[gi],
+ GaussPoints[gi],
+ magSi,
+ magSj,
+ di,
+ dj,
+ ci,
+ cj,
+ cosij,
+ alpha,
+ quadOrder
+ );
+
+ F2LIij += dIntFij;
+ }
+ }
+
+ scalar err = (F2LIij-oldEToeInt)/F2LIij;
+
+ if
+ (
+ (mag(err) < tolGaussQuad && gi > 0)
+ || gi == maxQuadOrder-1
+ )
+ {
+ F2LI[coarseFaceI][visCoarseFaceI] =
+ F2LIij/mag(Ai)/4/constant::mathematical::pi;
+ break;
+ }
+ else
+ {
+ oldEToeInt = F2LIij;
+ }
+ }
+ }
- F[coarseFaceI][visCoarseFaceI] = Fij;
- sumViewFactorPatch[fromPatchId][toPatchId] += Fij*mag(Ai);
+ sumViewFactorPatch[fromPatchId][toPatchId] +=
+ F2LI[coarseFaceI][visCoarseFaceI]*mag(Ai);
}
}
}
-
- if (Pstream::master())
+ else
{
- Info << "Writing view factor matrix..." << endl;
+ FatalErrorInFunction
+ << " View factors are not available in 2D "
+ << exit(FatalError);
}
// Write view factors matrix in listlist form
- F.write();
+ F2LI.write();
reduce(sumViewFactorPatch, sumOp<scalarSquareMatrix>());
reduce(patchArea, sumOp<scalarList>());
-
if (Pstream::master() && debug)
{
forAll(viewFactorsPatches, i)
{
label patchI = viewFactorsPatches[i];
- forAll(viewFactorsPatches, i)
+ for (label j=i; j<viewFactorsPatches.size(); j++)
{
- label patchJ = viewFactorsPatches[i];
+ label patchJ = viewFactorsPatches[j];
+
Info << "F" << patchI << patchJ << ": "
<< sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
<< endl;
@@ -758,9 +1118,13 @@ int main(int argc, char *argv[])
}
}
-
if (writeViewFactors)
{
+ if (Pstream::master())
+ {
+ Info << "Writing view factor matrix..." << endl;
+ }
+
volScalarField viewFactorField
(
IOobject
@@ -791,13 +1155,14 @@ int main(int argc, char *argv[])
forAll(coarseToFine, coarseI)
{
- const scalar Fij = sum(F[compactI]);
+ scalar FiSum = sum(F2LI[compactI]);
+
const label coarseFaceID = coarsePatchFace[coarseI];
const labelList& fineFaces = coarseToFine[coarseFaceID];
forAll(fineFaces, fineId)
{
const label faceID = fineFaces[fineId];
- vfbf[patchID][faceID] = Fij;
+ vfbf[patchID][faceID] = FiSum;
}
compactI++;
}
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/Make/files b/applications/utilities/preProcessing/viewFactorsGenExt/Make/files
deleted file mode 100644
index 6d18b01b2c5207a925b59646e814faa9acfc122f..0000000000000000000000000000000000000000
--- a/applications/utilities/preProcessing/viewFactorsGenExt/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-viewFactorsGenExt.C
-
-EXE = $(FOAM_APPBIN)/viewFactorsGenExt
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/Make/options b/applications/utilities/preProcessing/viewFactorsGenExt/Make/options
deleted file mode 100644
index 2c9476be08fb6208ae96cf7f7d3af6972dea4bc1..0000000000000000000000000000000000000000
--- a/applications/utilities/preProcessing/viewFactorsGenExt/Make/options
+++ /dev/null
@@ -1,39 +0,0 @@
-PBRT = /home/sergio/pub/pbrt
-PBRTSRC = $(PBRT)/pbrt-v3/src
-EXTSRC = $(PBRT)/pbrt-v3/src/ext
-GLOGSRC = $(EXTSRC)/glog/src
-
-PBRTBUILD = /home/sergio/pub/pbrt/pbrt-v3-build
-EXTBUILD = $(PBRTBUILD)/src/ext
-GLOGBUILD = $(EXTBUILD)/glog
-
-EXE_INC = \
- -Wno-old-style-cast \
- -I$(PBRTSRC)/core \
- -I$(PBRTSRC)/accelerators \
- -I$(PBRTSRC)/shapes \
- -I$(PBRTSRC) \
- -I$(EXTSRC) \
- -I$(GLOGSRC) \
- -I$(GLOGBUILD) \
- -I./rays \
- -I$(LIB_SRC)/finiteVolume/lnInclude \
- -I$(LIB_SRC)/surfMesh/lnInclude \
- -I$(LIB_SRC)/meshTools/lnInclude \
- -I$(LIB_SRC)/parallel/distributed/lnInclude \
-
-EXE_LIBS = \
- -L$(GLOGBUILD) -lglog \
- -L$(PBRTBUILD) -lpbrt \
- -L$(EXTBUILD)/openexr/OpenEXR/IlmImf -lIlmImf \
- -L$(EXTBUILD)/openexr/IlmBase/Imath -lImath \
- -L$(EXTBUILD)/openexr/IlmBase/Half -lHalf \
- -L$(EXTBUILD)/openexr/IlmBase/Iex -lIex \
- -L$(EXTBUILD)/openexr/IlmBase/IexMath -lIexMath \
- -L$(EXTBUILD)/openexr/IlmBase/IlmThread -lIlmThread \
- -L$(EXTBUILD)/ptex/src/ptex -lPtex \
- -lfiniteVolume \
- -lsurfMesh \
- -lmeshTools \
- -ldistributed \
- -lradiationModels
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/searchingEngineExt.H b/applications/utilities/preProcessing/viewFactorsGenExt/searchingEngineExt.H
deleted file mode 100644
index 678899fc7dc182ff05ea17868f79b8558c780a66..0000000000000000000000000000000000000000
--- a/applications/utilities/preProcessing/viewFactorsGenExt/searchingEngineExt.H
+++ /dev/null
@@ -1,121 +0,0 @@
-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);
-
-labelList triSurfaceToAgglom(5*nFineFaces);
-
-const triSurface localSurface = triangulate
-(
- patches,
- includePatches,
- finalAgglom,
- triSurfaceToAgglom,
- globalNumbering,
- coarsePatches
-);
-
-
-distributedTriSurfaceMesh surfacesMesh
-(
- IOobject
- (
- "wallSurface.stl",
- runTime.constant(), // directory
- "triSurface", // instance
- runTime, // registry
- IOobject::NO_READ,
- IOobject::AUTO_WRITE
- ),
- localSurface,
- dict
-);
-
-
-triSurfaceToAgglom.resize(surfacesMesh.size());
-
-// pbrt surface
-std::vector<Point3f> vertices;
-
-const pointField& pts = surfacesMesh.localPoints();
-for (const auto& pt : pts)
-{
- vertices.push_back(Point3f(pt[0], pt[1], pt[2]));
-}
-
-std::vector<int> indices;
-for (const auto& tri : surfacesMesh.localFaces())
-{
- indices.push_back(tri[0]);
- indices.push_back(tri[1]);
- indices.push_back(tri[2]);
-}
-
-std::vector<int> faceIndices;
-faceIndices.reserve(surfacesMesh.localFaces().size());
-
-for (label faceI = 0; faceI < surfacesMesh.localFaces().size(); faceI++)
-{
- faceIndices.push_back(faceI);
-}
-
-Transform o2w;
-Transform w2o;
-
-std::vector<std::shared_ptr<Shape>> surfacesMesh_pbrt
-(
- CreateTriangleMesh
- (
- &o2w,
- &w2o,
- false, // bool reverseOrientation
- surfacesMesh.size(), // int nTriangles
- &indices[0], // int *vertexIndices
- vertices.size(),
- &vertices[0],
- nullptr,
- nullptr,
- nullptr,
- nullptr,
- nullptr,
- &faceIndices[0]
- )
-);
-
-std::vector<std::shared_ptr<Primitive>> prims;
-
-prims.reserve(surfacesMesh_pbrt.size());
-for (auto s : surfacesMesh_pbrt)
-{
- prims.push_back
- (
- std::make_shared<GeometricPrimitive>
- (s, nullptr, nullptr, nullptr)
- );
-}
-
-std::shared_ptr<Primitive> accel;
-
-ParamSet paramSet;
-
-accel = CreateBVHAccelerator(std::move(prims), paramSet);
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H b/applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H
deleted file mode 100644
index 0cc8240d9354bcced0d181863f2d1d015ba727cc..0000000000000000000000000000000000000000
--- a/applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H
+++ /dev/null
@@ -1,196 +0,0 @@
-// All rays expressed as start face (local) index end end face (global)
-// Pre-size by assuming a certain percentage is visible.
-
-// Maximum length for dynamicList
-const label maxDynListLength
-(
- viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000000)
-);
-
-for (const int proci : Pstream::allProcs())
-{
-
- std::vector<pbrt::Point3f> start;
- start.reserve(coarseMesh.nFaces());
- std::vector<pbrt::Vector3f> dir;
- dir.reserve(start.size());
-
- std::vector<pbrt::Point3f> end(start.size());
- end.reserve(start.size());
-
- DynamicList<label> startIndex(start.size());
- DynamicList<label> endIndex(start.size());
-
- DynamicList<label> startAgg(start.size());
- DynamicList<label> endAgg(start.size());
-
- const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
- const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
- const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
-
- const pointField& remoteArea = remoteCoarseSf[proci];
- const pointField& remoteFc = remoteCoarseCf[proci];
- const labelField& remoteAgg = remoteCoarseAgg[proci];
-
- pbrt::Vector3f smallV(SMALL, SMALL, SMALL);
-
- label i = 0;
- label j = 0;
- do
- {
- for (; i < myFc.size(); i++)
- {
- const point& fc = myFc[i];
- const vector& fA = myArea[i];
- const label& fAgg = myAgg[i];
-
- for (; j < remoteFc.size(); j++)
- {
- if (proci != Pstream::myProcNo() || i != j)
- {
- const point& remFc = remoteFc[j];
- const vector& remA = remoteArea[j];
- const label& remAgg = remoteAgg[j];
-
- const vector d = remFc - fc;
-
- const vector nd = d/mag(d);
- const vector nfA = fA/mag(fA);
- const vector nremA = remA/mag(remA);
-
- if (((nd & nfA) < 0) && ((nd & nremA) > 0))
- {
- pbrt::Vector3f direction(d[0], d[1], d[2]);
- direction += smallV;
- pbrt::Point3f s(fc[0], fc[1], fc[2]);
- end.push_back(s + direction);
-
- s += pbrt::Point3f(0.01*d[0], 0.01*d[1], 0.01*d[2]);
-
- start.push_back(s);
- startIndex.append(i);
- startAgg.append(globalNumbering.toGlobal(proci, fAgg));
-
- dir.push_back(direction);
-
- label globalI = globalNumbering.toGlobal(proci, j);
- endIndex.append(globalI);
- endAgg.append(globalNumbering.toGlobal(proci, remAgg));
- if (startIndex.size() > maxDynListLength)
- {
- FatalErrorInFunction
- << "Dynamic list need from capacity."
- << "Actual size maxDynListLength : "
- << maxDynListLength
- << abort(FatalError);
- }
- }
- }
- }
-
- if (j == remoteFc.size())
- {
- j = 0;
- }
- }
-
- } while (i < myFc.size());
-
- label totalRays(startIndex.size());
- reduce(totalRays, sumOp<scalar>());
-
- if (Pstream::master() && debug)
- {
- Info<< "Number of rays :" << totalRays << endl;
- }
-
- DynamicList<label> dRayIs;
-
- std::vector<pbrt::Ray> raysInt;
- raysInt.reserve(start.size());
-
- std::vector<int> raysTriIndex;
- raysTriIndex.reserve(start.size());
-
- for (unsigned long rayI = 0; rayI < start.size(); ++rayI)
- {
- pbrt::Ray ray(start[rayI], dir[rayI]);
-
- pbrt::SurfaceInteraction isect;
- bool intersects = accel->Intersect(ray, &isect);
- const Triangle* tri = dynamic_cast<const Triangle *>(isect.shape);
-
- if (intersects)
- {
- const int index = tri->faceIndex;
-
- const Vector3f delta(ray(ray.tMax) - end[rayI]);
-
- if (delta.Length() < intTol)
- {
- rayStartFace.append(startIndex[rayI]);
- rayEndFace.append(endIndex[rayI]);
- }
- else if (triSurfaceToAgglom[index] == startAgg[rayI])
- {
- dRayIs.append(rayI);
- raysInt.push_back(ray);
- raysTriIndex.push_back(index);
- }
- }
- }
- start.clear();
-
- labelField rayIs;
- rayIs.transfer(dRayIs);
- dRayIs.clear();
-
- boolList converged(rayIs.size(), false);
- label nConverged = 0;
- label iter = 0;
- do
- {
- forAll(rayIs, rayI)
- {
- const label rayID = rayIs[rayI];
- pbrt::Ray& rayHit = raysInt[rayI];
- const int index = raysTriIndex[rayI];
-
- if (!converged[rayI])
- {
- if (triSurfaceToAgglom[index] == startAgg[rayID])
- {
- rayHit.tMax *= 1.1;
- rayHit.o = rayHit(rayHit.tMax);
-
- pbrt::SurfaceInteraction isect;
- bool intersects = accel->Intersect(rayHit, &isect);
- if (intersects)
- {
- const Triangle* tri =
- dynamic_cast<const Triangle *>(isect.shape);
- const int index = tri->faceIndex;
-
- raysTriIndex[rayI] = index;
- }
- }
- else if (triSurfaceToAgglom[index] == endAgg[rayID])
- {
- converged[rayI] = true;
- nConverged++;
- rayStartFace.append(startIndex[rayID]);
- rayEndFace.append(endIndex[rayID]);
- }
- }
- }
-
- iter ++;
-
- } while (nConverged < converged.size() && iter < 10);
-
- startIndex.clear();
- end.clear();
- endIndex.clear();
- startAgg.clear();
- endAgg.clear();
-}
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C b/applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C
deleted file mode 100644
index 2f5f9ba496c089201072330b2dbbf003eb04880d..0000000000000000000000000000000000000000
--- a/applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C
+++ /dev/null
@@ -1,1227 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | www.openfoam.com
- \\/ M anipulation |
--------------------------------------------------------------------------------
- Copyright (C) 2022 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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
- viewFactorsGenExt
-
-Group
- grpPreProcessingUtilities
-
-Description
- This view factors generation application uses a combined approach of
- double area integral (2AI) and double linear integral (2LI). 2AI is used
- when the two surfaces are 'far' apart and 2LI whenre they are 'close'.
- 2LI is integrated along edges using Gaussian quadrature.
- The distance between faces is calculating a ratio between averaged areas
- and the distance between face centres.
-
- The input from viewFactorsDict are:
-
- tolGaussQuad 0.1; // GaussQuad error
- distTol 8; // R/Average(rm)
- alpha 0.22; // Use for common edges for 2LI
-
-
- For debugging purposes, the following entries can be set in viewFactorsDict:
-
- writeViewFactorMatrix true;
- writeFacesAgglomeration false;
- dumpRays false;
-
-
- writeViewFactorMatrix writes the sum of the VF on each face.
- writeFacesAgglomeration writes the agglomeration
- dumpRays dumps rays
-
-
- In order to specify the participants patches in the VF calculation the
- keywaord viewFactorWall should be added to the boundary file.
-
- floor
- {
- type wall;
- inGroups 2(wall viewFactorWall);
- nFaces 100;
- startFace 3100;
- }
-
-\*---------------------------------------------------------------------------*/
-
-#include "argList.H"
-#include "fvMesh.H"
-#include "volFields.H"
-#include "surfaceFields.H"
-#include "distributedTriSurfaceMesh.H"
-#include "meshTools.H"
-#include "constants.H"
-
-#include "uindirectPrimitivePatch.H"
-#include "DynamicField.H"
-#include "unitConversion.H"
-
-#include "scalarMatrices.H"
-#include "labelListIOList.H"
-#include "scalarListIOList.H"
-
-#include "singleCellFvMesh.H"
-
-#include "IOmapDistribute.H"
-
-#define PBRT_CONSTEXPR constexpr
-#define PBRT_THREAD_LOCAL thread_local
-
-#include "pbrt.h"
-#include "shape.h"
-#include "triangle.h"
-#include "geometry.h"
-#include "paramset.h"
-
-#include "accelerators/bvh.h"
-
-using namespace Foam;
-using namespace Foam::constant;
-using namespace Foam::constant::mathematical;
-using namespace pbrt;
-
-
-triSurface triangulate
-(
- const polyBoundaryMesh& bMesh,
- const labelHashSet& includePatches,
- const labelListIOList& finalAgglom,
- labelList& triSurfaceToAgglom,
- const globalIndex& globalNumbering,
- const polyBoundaryMesh& coarsePatches
-)
-{
- const polyMesh& mesh = bMesh.mesh();
-
- // Storage for surfaceMesh. Size estimate.
- DynamicList<labelledTri> triangles(mesh.nBoundaryFaces());
-
- label newPatchI = 0;
- label localTriFaceI = 0;
-
- for (const label patchI : includePatches)
- {
- const polyPatch& patch = bMesh[patchI];
- const pointField& points = patch.points();
-
- label nTriTotal = 0;
-
- forAll(patch, patchFaceI)
- {
- const face& f = patch[patchFaceI];
-
- faceList triFaces(f.nTriangles(points));
-
- label nTri = 0;
-
- f.triangles(points, nTri, triFaces);
-
- forAll(triFaces, triFaceI)
- {
- const face& f = triFaces[triFaceI];
-
- triangles.append(labelledTri(f[0], f[1], f[2], newPatchI));
-
- nTriTotal++;
-
- triSurfaceToAgglom[localTriFaceI++] = globalNumbering.toGlobal
- (
- Pstream::myProcNo(),
- finalAgglom[patchI][patchFaceI]
- + coarsePatches[patchI].start()
- );
- }
- }
-
- newPatchI++;
- }
-
- //triSurfaceToAgglom.resize(localTriFaceI-1);
-
- triangles.shrink();
-
- triSurface rawSurface(triangles, mesh.points());
-
- triSurface surface
- (
- rawSurface.localFaces(),
- rawSurface.localPoints()
- );
-
- // Combine the triSurfaces across all processors
- if (Pstream::parRun())
- {
- List<List<labelledTri>> surfaceProcTris(Pstream::nProcs());
- List<pointField> surfaceProcPoints(Pstream::nProcs());
-
- surfaceProcTris[Pstream::myProcNo()] = surface;
- surfaceProcPoints[Pstream::myProcNo()] = surface.points();
-
- Pstream::gatherList(surfaceProcTris);
- Pstream::scatterList(surfaceProcTris);
- Pstream::gatherList(surfaceProcPoints);
- Pstream::scatterList(surfaceProcPoints);
-
- label nTris = 0;
- forAll(surfaceProcTris, i)
- {
- nTris += surfaceProcTris[i].size();
- }
-
- List<labelledTri> globalSurfaceTris(nTris);
- label trii = 0;
- label offset = 0;
- forAll(surfaceProcTris, i)
- {
- forAll(surfaceProcTris[i], j)
- {
- globalSurfaceTris[trii] = surfaceProcTris[i][j];
- globalSurfaceTris[trii][0] += offset;
- globalSurfaceTris[trii][1] += offset;
- globalSurfaceTris[trii][2] += offset;
- trii++;
- }
- offset += surfaceProcPoints[i].size();
- }
-
- label nPoints = 0;
- forAll(surfaceProcPoints, i)
- {
- nPoints += surfaceProcPoints[i].size();
- }
-
- pointField globalSurfacePoints(nPoints);
-
- label pointi = 0;
- forAll(surfaceProcPoints, i)
- {
- forAll(surfaceProcPoints[i], j)
- {
- globalSurfacePoints[pointi++] = surfaceProcPoints[i][j];
- }
- }
-
- surface = triSurface(globalSurfaceTris, globalSurfacePoints);
- }
-
- // Add patch names to surface
- surface.patches().setSize(newPatchI);
-
- newPatchI = 0;
-
- for (const label patchI : includePatches)
- {
- const polyPatch& patch = bMesh[patchI];
-
- surface.patches()[newPatchI].index() = patchI;
- surface.patches()[newPatchI].name() = patch.name();
- surface.patches()[newPatchI].geometricType() = patch.type();
-
- newPatchI++;
- }
-
- return surface;
-}
-
-
-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;
- }
- }
- str.flush();
-}
-
-
-scalar calculateViewFactorFij2AI
-(
- const vector& i,
- const vector& j,
- const vector& dAi,
- const vector& dAj
-)
-{
- vector r = i - j;
- scalar rMag = mag(r);
-
- if (rMag > SMALL)
- {
- 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)
- );
- }
- else
- {
- return 0;
- }
-}
-
-
-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];
- }
- }
-}
-
-scalar GaussQuad
-(
-
- const scalarList& w,
- const scalarList& p,
- const scalar& magSi,
- const scalar& magSj,
- const vector& di,
- const vector& dj,
- const vector& ci,
- const vector& cj,
- const scalar cosij,
- const scalar alpha,
- label gi
-)
-{
- scalar dIntFij = 0;
- if (gi == 0)
- {
- vector r(ci - cj);
- if (mag(r) < SMALL)
- {
- r = (alpha*magSi)*di;
- }
- dIntFij = max(cosij*Foam::log(r&r)*magSi*magSj, 0);
- }
- else
- {
- List<vector> pi(w.size());
- forAll (pi, i)
- {
- pi[i] = ci + p[i]*(magSi/2)*di;
- }
-
-
- List<vector> pj(w.size());
- forAll (pj, i)
- {
- pj[i] = cj + p[i]*(magSj/2)*dj;
- }
-
- forAll (w, i)
- {
- forAll (w, j)
- {
- vector r(pi[i] - pj[j]);
-
- if (mag(r) < SMALL)
- {
- r = (alpha*magSi)*di;
- dIntFij +=
- cosij*w[i]*w[j]*Foam::log(r&r);
- }
- else
- {
- dIntFij +=
- cosij*w[i]*w[j]*Foam::log(r&r);
- }
-
- }
- }
-
- dIntFij *= (magSi/2) * (magSj/2);
-
- }
- return dIntFij;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
- argList::addNote
- (
- "Calculate view factors from face agglomeration array."
- " The finalAgglom generated by faceAgglomerate utility."
- );
-
- #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 word viewFactorWall("viewFactorWall");
-
- const bool writeViewFactors =
- viewFactorDict.getOrDefault("writeViewFactorMatrix", false);
-
- const bool dumpRays =
- viewFactorDict.getOrDefault("dumpRays", false);
-
- const label debug = viewFactorDict.getOrDefault<label>("debug", 0);
-
- const scalar tolGaussQuad =
- viewFactorDict.getOrDefault<scalar>("tolGaussQuad", 0.01);
-
- const scalar distTol =
- viewFactorDict.getOrDefault<scalar>("distTol", 8);
-
- const scalar alpha =
- viewFactorDict.getOrDefault<scalar>("alpha", 0.21);
-
- const scalar intTol =
- viewFactorDict.getOrDefault<scalar>("intTol", 1e-2);
-
- // 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)
- {
- Pout << "\nCreating single cell mesh..." << endl;
- }
-
- singleCellFvMesh coarseMesh
- (
- IOobject
- (
- "coarse:" + mesh.name(),
- runTime.timeName(),
- runTime,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- mesh,
- finalAgglom
- );
-
- if (debug)
- {
- Pout << "\nCreated single cell mesh..." << endl;
- }
-
-
- // 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.indices(viewFactorWall));
- for (const label patchi : viewFactorsPatches)
- {
- nCoarseFaces += coarsePatches[patchi].size();
- nFineFaces += patches[patchi].size();
- }
-
- // 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);
- DynamicList<label> localAgg(nCoarseFaces);
- labelHashSet includePatches;
-
- for (const label patchID : viewFactorsPatches)
- {
- const polyPatch& pp = patches[patchID];
- const labelList& agglom = finalAgglom[patchID];
-
- includePatches.insert(patchID);
-
- if (agglom.size() > 0)
- {
- 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];
- const label agglomI =
- agglom[fineFaces[0]] + coarsePatches[patchID].start();
-
- // 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()
- ) = upp.faceCentres();
-
- SubList<point>
- (
- availablePoints,
- upp.localPoints().size(),
- upp.faceCentres().size()
- ) = 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);
- localAgg.append(agglomI);
-
- }
- }
- }
-
- // Distribute local coarse Cf and Sf for shooting rays
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- List<pointField> remoteCoarseCf(Pstream::nProcs());
- List<pointField> remoteCoarseSf(Pstream::nProcs());
- List<labelField> remoteCoarseAgg(Pstream::nProcs());
-
- remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
- remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
- remoteCoarseAgg[Pstream::myProcNo()] = localAgg;
-
- Pstream::gatherList(remoteCoarseCf);
- Pstream::scatterList(remoteCoarseCf);
- Pstream::gatherList(remoteCoarseSf);
- Pstream::scatterList(remoteCoarseSf);
- Pstream::gatherList(remoteCoarseAgg);
- Pstream::scatterList(remoteCoarseAgg);
-
-
- globalIndex globalNumbering(nCoarseFaces);
-
- // Set up searching engine for obstacles
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- #include "searchingEngineExt.H"
-
- // Determine rays between coarse face centres
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- DynamicList<label> rayStartFace(nCoarseFaces + 0.01*nCoarseFaces);
-
- DynamicList<label> rayEndFace(rayStartFace.size());
-
- // Return rayStartFace in local index and rayEndFace in global index
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- #include "shootRaysExt.H"
-
- // Calculate number of visible faces from local index
- labelList nVisibleFaceFaces(nCoarseFaces, Zero);
-
- 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);
-
- // 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
- // (2AA) need coarse (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
- // (2LI) need edges (li)
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- pointField compactCoarseCf(map.constructSize(), Zero);
- pointField compactCoarseSf(map.constructSize(), Zero);
- List<List<point>> compactFineSf(map.constructSize());
- List<List<point>> compactFineCf(map.constructSize());
-
- DynamicList<List<point>> compactPoints(map.constructSize());
-
- DynamicList<label> compactPatchId(map.constructSize());
-
- // Insert my coarse local values
- SubList<point>(compactCoarseSf, nCoarseFaces) = localCoarseSf;
- SubList<point>(compactCoarseCf, nCoarseFaces) = localCoarseCf;
-
- const faceList& faces = mesh.faces();
-
- // Insert my fine local values
- label compactI = 0;
- forAll(viewFactorsPatches, i)
- {
- label patchID = viewFactorsPatches[i];
-
- const labelList& agglom = finalAgglom[patchID];
- if (agglom.size() > 0)
- {
- label nAgglom = max(agglom)+1;
- labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
- const labelList& coarsePatchFace =
- coarseMesh.patchFaceMap()[patchID];
-
- const polyPatch& pp = patches[patchID];
-
- forAll(coarseToFine, coarseI)
- {
- compactPatchId.append(patchID);
- List<point>& fineCf = compactFineCf[compactI];
- List<point>& fineSf = compactFineSf[compactI];
-
- label startFace = pp.start();
-
- const vectorField locPoints
- (
- mesh.points(),
- faces[coarseI + startFace]
- );
-
- const label coarseFaceI = coarsePatchFace[coarseI];
- const labelList& fineFaces = coarseToFine[coarseFaceI];
-
- fineCf.setSize(fineFaces.size());
- fineSf.setSize(fineFaces.size());
-
- compactPoints.append(locPoints);
-
- fineCf = UIndirectList<point>
- (
- mesh.Cf().boundaryField()[patchID],
- coarseToFine[coarseFaceI]
- );
- fineSf = UIndirectList<point>
- (
- mesh.Sf().boundaryField()[patchID],
- coarseToFine[coarseFaceI]
- );
-
- compactI++;
- }
- }
- }
-
- if (Pstream::master())
- {
- Info<< "map distribute..." << endl;
- }
-
- // Do all swapping
- map.distribute(compactCoarseSf);
- map.distribute(compactCoarseCf);
- map.distribute(compactFineCf);
- map.distribute(compactFineSf);
-
- map.distribute(compactPoints);
-
- 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 F2LI
- (
- 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,
- 0.0
- );
-
- scalarList patchArea(totalPatches, Zero);
-
- if (Pstream::master())
- {
- Info<< "\nCalculating view factors..." << endl;
- }
-
- List<scalarList> GaussPoints(5);
- GaussPoints[0].setSize(1);
- GaussPoints[0] = 0;
-
- GaussPoints[1].setSize(2);
- GaussPoints[1][0] = 1/std::sqrt(3);
- GaussPoints[1][1] = -1/std::sqrt(3);
-
- GaussPoints[2].setSize(3);
- GaussPoints[2][0] = 0;
- GaussPoints[2][1] = 0.774597;
- GaussPoints[2][2] = -0.774597;
-
- GaussPoints[3].setSize(4);
- GaussPoints[3][0] = 0.339981;
- GaussPoints[3][1] = -0.339981;
- GaussPoints[3][2] = 0.861136;
- GaussPoints[3][3] = -0.861136;
-
- GaussPoints[4].setSize(5);
- GaussPoints[4][0] = 0;
- GaussPoints[4][1] = 0.538469;
- GaussPoints[4][2] = -0.538469;
- GaussPoints[4][3] = 0.90618;
- GaussPoints[4][4] = -0.90618;
-
-
- List<scalarList> GaussWeights(5);
- GaussWeights[0].setSize(1);
- GaussWeights[0] = 2;
-
- GaussWeights[1].setSize(2);
- GaussWeights[1][0] = 1;
- GaussWeights[1][1] = 1;
-
- GaussWeights[2].setSize(3);
- GaussWeights[2][0] = 0.888889;
- GaussWeights[2][1] = 0.555556;
- GaussWeights[2][2] = 0.555556;
-
- GaussWeights[3].setSize(4);
- GaussWeights[3][0] = 0.652145;
- GaussWeights[3][1] = 0.652145;
- GaussWeights[3][2] = 0.347855;
- GaussWeights[3][3] = 0.347855;
-
- GaussWeights[4].setSize(5);
- GaussWeights[4][0] = 0.568889;
- GaussWeights[4][1] = 0.478629;
- GaussWeights[4][2] = 0.478629;
- GaussWeights[4][3] = 0.236927;
- GaussWeights[4][4] = 0.236927;
-
- label maxQuadOrder = 5;
-
- 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];
-
- const List<point>& lPoints = compactPoints[coarseFaceI];
-
- patchArea[fromPatchId] += mag(Ai);
-
- const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
-
- forAll(visCoarseFaces, visCoarseFaceI)
- {
- //F2AI[coarseFaceI].setSize(visCoarseFaces.size());
- F2LI[coarseFaceI].setSize(visCoarseFaces.size());
- label compactJ = visCoarseFaces[visCoarseFaceI];
- const List<point>& remoteFineSj = compactFineSf[compactJ];
- const List<point>& remoteFineCj = compactFineCf[compactJ];
-
- const List<point>& rPointsCj = compactPoints[compactJ];
-
- const label toPatchId = compactPatchId[compactJ];
-
- bool far(false);
- // Relative distance
- forAll(localFineSf, i)
- {
- const scalar dAi =
- Foam::sqrt
- (
- mag(localFineSf[i])/constant::mathematical::pi
- );
- const vector& dCi = localFineCf[i];
-
- forAll(remoteFineSj, j)
- {
- const scalar dAj =
- Foam::sqrt
- (
- mag(remoteFineSj[j])/constant::mathematical::pi
- );
- const vector& dCj = remoteFineCj[j];
-
- const scalar dist = mag(dCi - dCj)/((dAi + dAj)/2);
-
- if (dist > distTol)
- {
- far = true;
- }
- }
- }
-
- if (far)
- {
- // 2AI method
- scalar F2AIij = 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 = calculateViewFactorFij2AI
- (
- dCi,
- dCj,
- dAi,
- dAj
- );
-
- F2AIij += dIntFij;
- }
- }
- F2LI[coarseFaceI][visCoarseFaceI] = F2AIij/mag(Ai);
- }
- else
- {
- // 2LI method
- label nLocal = lPoints.size();
- label nRemote = rPointsCj.size();
-
- // Using sub-divisions (quadrature)
- scalar oldEToeInt = 0;
- for (label gi=0; gi < maxQuadOrder; gi++)
- {
- scalar F2LIij = 0;
- for(label i=0; i<nLocal; i++)
- {
- vector si;
- vector ci;
-
- vector sj;
- vector cj;
-
- if (i == 0)
- {
- si = lPoints[i] - lPoints[nLocal-1];
- ci = (lPoints[i] + lPoints[nLocal-1])/2;
- }
- else
- {
- si = lPoints[i] - lPoints[i-1];
- ci = (lPoints[i] + lPoints[i-1])/2;
- }
-
- for(label j=0; j<nRemote; j++)
- {
- if (j == 0)
- {
- sj = rPointsCj[j]-rPointsCj[nRemote-1];
- cj = (rPointsCj[j]+rPointsCj[nRemote-1])/2;
- }
- else
- {
- sj = rPointsCj[j] - rPointsCj[j-1];
- cj = (rPointsCj[j] + rPointsCj[j-1])/2;
- }
-
-
- scalar magSi = mag(si);
- scalar magSj = mag(sj);
- scalar cosij = (si & sj)/(magSi * magSj);
-
- vector di = si/magSi;
- vector dj = sj/magSj;
-
- label quadOrder = gi;
- const vector r(ci - cj);
- // Common edges use n = 0
- if (mag(r) < SMALL)
- {
- quadOrder = 0;
- }
-
- scalar dIntFij =
- GaussQuad
- (
- GaussWeights[gi],
- GaussPoints[gi],
- magSi,
- magSj,
- di,
- dj,
- ci,
- cj,
- cosij,
- alpha,
- quadOrder
- );
-
- F2LIij += dIntFij;
- }
- }
-
- scalar err = (F2LIij-oldEToeInt)/F2LIij;
-
- if
- (
- (mag(err) < tolGaussQuad && gi > 0)
- || gi == maxQuadOrder-1
- )
- {
- //DebugVar(gi)
- F2LI[coarseFaceI][visCoarseFaceI] =
- F2LIij/mag(Ai)/4/constant::mathematical::pi;
- break;
- }
- else
- {
- oldEToeInt = F2LIij;
- }
- }
- }
-
- sumViewFactorPatch[fromPatchId][toPatchId] +=
- F2LI[coarseFaceI][visCoarseFaceI]*mag(Ai);
- }
- }
- }
- else
- {
- FatalErrorInFunction
- << " View factors are not available in 2D "
- << exit(FatalError);
- }
-
-
-
- // Write view factors matrix in listlist form
- F2LI.write();
-
- reduce(sumViewFactorPatch, sumOp<scalarSquareMatrix>());
- reduce(patchArea, sumOp<scalarList>());
-
- if (Pstream::master() && debug)
- {
- forAll(viewFactorsPatches, i)
- {
- label patchI = viewFactorsPatches[i];
- for (label j=i; j<viewFactorsPatches.size(); j++)
- {
- label patchJ = viewFactorsPatches[j];
-
- Info << "F" << patchI << patchJ << ": "
- << sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
- << endl;
- }
- }
- }
-
- if (writeViewFactors)
- {
- if (Pstream::master())
- {
- Info << "Writing view factor matrix..." << endl;
- }
-
- volScalarField viewFactorField
- (
- IOobject
- (
- "viewFactorField",
- mesh.time().timeName(),
- mesh,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- mesh,
- dimensionedScalar(dimless, Zero)
- );
-
- label compactI = 0;
-
- volScalarField::Boundary& vfbf = viewFactorField.boundaryFieldRef();
- forAll(viewFactorsPatches, i)
- {
- label patchID = viewFactorsPatches[i];
- const labelList& agglom = finalAgglom[patchID];
- if (agglom.size() > 0)
- {
- label nAgglom = max(agglom)+1;
- labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
- const labelList& coarsePatchFace =
- coarseMesh.patchFaceMap()[patchID];
-
- forAll(coarseToFine, coarseI)
- {
- scalar FiSum = sum(F2LI[compactI]);
-
- const label coarseFaceID = coarsePatchFace[coarseI];
- const labelList& fineFaces = coarseToFine[coarseFaceID];
- forAll(fineFaces, fineId)
- {
- const label faceID = fineFaces[fineId];
- vfbf[patchID][faceID] = FiSum;
- }
- 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];
-
- forAllConstIters(localToCompactMap, iter)
- {
- compactToGlobal[*iter] = globalNumbering.toGlobal
- (
- procI,
- iter.key()
- );
- }
- }
-
-
- 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
- ),
- std::move(globalFaceFaces)
- );
- IOglobalFaceFaces.write();
-
-
- IOmapDistribute IOmapDist
- (
- IOobject
- (
- "mapDist",
- mesh.facesInstance(),
- mesh,
- IOobject::NO_READ,
- IOobject::AUTO_WRITE
- ),
- std::move(map)
- );
-
- IOmapDist.write();
-
- Info<< "End\n" << endl;
- return 0;
-}
-
-
-// ************************************************************************* //