From c08fc5ecd9901dfb1df1fddc28a2107925c34f37 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 4 Aug 2022 14:24:50 +0100
Subject: [PATCH] ENH: viewFactorGen: revert to v2206 without CGAL
---
.../preProcessing/viewFactorsGen/Allwmake | 19 +++
.../preProcessing/viewFactorsGen/Make/options | 3 +-
.../viewFactorsGen/searchingEngine.H | 31 +---
.../viewFactorsGen/searchingEngine_CGAL.H | 84 ++++++++++
.../preProcessing/viewFactorsGen/shootRays.H | 145 +++++++++++++-----
.../viewFactorsGen/shootRays_CGAL.H | 131 ++++++++++++++++
.../viewFactorsGen/viewFactorsGen.C | 41 +++--
7 files changed, 370 insertions(+), 84 deletions(-)
create mode 100755 applications/utilities/preProcessing/viewFactorsGen/Allwmake
create mode 100644 applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H
create mode 100644 applications/utilities/preProcessing/viewFactorsGen/shootRays_CGAL.H
diff --git a/applications/utilities/preProcessing/viewFactorsGen/Allwmake b/applications/utilities/preProcessing/viewFactorsGen/Allwmake
new file mode 100755
index 00000000000..f0d604d4e3e
--- /dev/null
+++ b/applications/utilities/preProcessing/viewFactorsGen/Allwmake
@@ -0,0 +1,19 @@
+#!/bin/sh
+cd "${0%/*}" || exit # Run from this directory
+. ${WM_PROJECT_DIR:?}/wmake/scripts/AllwmakeParseArguments # (error catching)
+. ${WM_PROJECT_DIR:?}/wmake/scripts/have_cgal
+
+#------------------------------------------------------------------------------
+unset COMP_FLAGS LINK_FLAGS
+
+if have_cgal
+then
+ echo " found CGAL -- enabling CGAL support."
+else
+ echo " did not find CGAL -- disabling CGAL functionality"
+ export COMP_FLAGS="-DNO_CGAL"
+fi
+
+wmake $targetType
+
+#------------------------------------------------------------------------------
diff --git a/applications/utilities/preProcessing/viewFactorsGen/Make/options b/applications/utilities/preProcessing/viewFactorsGen/Make/options
index 6ac3efad51f..08e64c8a7d4 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/Make/options
+++ b/applications/utilities/preProcessing/viewFactorsGen/Make/options
@@ -1,7 +1,8 @@
-include $(GENERAL_RULES)/cgal
+include $(GENERAL_RULES)/cgal-header-only
EXE_INC = \
-Wno-old-style-cast \
+ $(COMP_FLAGS) \
${CGAL_INC} \
-DCGAL_HEADER_ONLY \
-I$(LIB_SRC)/finiteVolume/lnInclude \
diff --git a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
index 667470d3f28..df38326ded0 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);
-triSurface localSurface = triangulate
+const triSurface localSurface = triangulate
(
patches,
includePatches,
@@ -36,8 +36,6 @@ triSurface localSurface = triangulate
);
-// CGAL surface
-
distributedTriSurfaceMesh surfacesMesh
(
IOobject
@@ -53,32 +51,7 @@ distributedTriSurfaceMesh surfacesMesh
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/searchingEngine_CGAL.H b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H
new file mode 100644
index 00000000000..667470d3f28
--- /dev/null
+++ b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H
@@ -0,0 +1,84 @@
+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);
+
+triSurface localSurface = triangulate
+(
+ patches,
+ includePatches,
+ finalAgglom,
+ triSurfaceToAgglom,
+ globalNumbering,
+ coarsePatches
+);
+
+
+// CGAL surface
+
+distributedTriSurfaceMesh surfacesMesh
+(
+ IOobject
+ (
+ "wallSurface.stl",
+ runTime.constant(), // directory
+ "triSurface", // instance
+ runTime, // registry
+ IOobject::NO_READ,
+ IOobject::NO_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 c738784b799..a5f808c3923 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
@@ -1,23 +1,35 @@
+// All rays expressed as start face (local) index end end face (global)
+// Pre-size by assuming a certain percentage is visible.
+
+
+if (Pstream::master())
+{
+ Info<< "\nShooting rays using distributedTriSurfaceMesh (no CGAL support)"
+ << endl;
+}
+
+
+
// Maximum length for dynamicList
const label maxDynListLength
(
- viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000000)
+ viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000)
);
for (const int proci : Pstream::allProcs())
{
- std::vector<Point> start;
- start.reserve(coarseMesh.nFaces());
-
- std::vector<Point> end(start.size());
- end.reserve(start.size());
+ // Shoot rays from me to proci. Note that even if processor has
+ // 0 faces we still need to call findLine to keep calls synced.
+ 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()];
@@ -46,31 +58,15 @@ for (const int proci : Pstream::allProcs())
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))
+ if (((d & fA) < 0.) && ((d & remA) > 0))
{
- 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]));
+ start.append(fc + 0.001*d);
startIndex.append(i);
- if (useAgglomeration)
- {
- startAgg.append(globalNumbering.toGlobal(proci, fAgg));
- endAgg.append(globalNumbering.toGlobal(proci, remAgg));
- }
-
+ startAgg.append(globalNumbering.toGlobal(proci, fAgg));
+ end.append(fc + 0.999*d);
label globalI = globalNumbering.toGlobal(proci, j);
endIndex.append(globalI);
-
+ endAgg.append(globalNumbering.toGlobal(proci, remAgg));
if (startIndex.size() > maxDynListLength)
{
FatalErrorInFunction
@@ -89,33 +85,98 @@ for (const int proci : Pstream::allProcs())
}
}
- } while (i < myFc.size());
+ } while (returnReduce(i < myFc.size(), orOp<bool>()));
- label totalRays(startIndex.size());
- reduce(totalRays, sumOp<scalar>());
+ List<pointIndexHit> hitInfo(startIndex.size());
+ surfacesMesh.findLine(start, end, hitInfo);
- if (debug)
- {
- Pout<< "Number of rays :" << totalRays << endl;
- }
+ // Return hit global agglo index
+ labelList aggHitIndex;
+ surfacesMesh.getField(hitInfo, aggHitIndex);
- for (unsigned long rayI = 0; rayI < start.size(); ++rayI)
- {
- Segment ray(start[rayI], end[rayI]);
- Segment_intersection intersects = tree.any_intersection(ray);
+ DynamicList<label> dRayIs;
- if (!intersects)
+ // 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)
+ {
+ if (!hitInfo[rayI].hit())
{
rayStartFace.append(startIndex[rayI]);
rayEndFace.append(endIndex[rayI]);
}
+ else if (aggHitIndex[rayI] == startAgg[rayI])
+ {
+ dRayIs.append(rayI);
+ }
}
+
start.clear();
- if (debug)
+
+ // 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
{
- Pout << "hits : " << rayStartFace.size()<< endl;
- }
+ 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);
startIndex.clear();
end.clear();
diff --git a/applications/utilities/preProcessing/viewFactorsGen/shootRays_CGAL.H b/applications/utilities/preProcessing/viewFactorsGen/shootRays_CGAL.H
new file mode 100644
index 00000000000..8be6e00c571
--- /dev/null
+++ b/applications/utilities/preProcessing/viewFactorsGen/shootRays_CGAL.H
@@ -0,0 +1,131 @@
+// Maximum length for dynamicList
+const label maxDynListLength
+(
+ viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000000)
+);
+
+for (const int proci : Pstream::allProcs())
+{
+ std::vector<Point> start;
+ start.reserve(coarseMesh.nFaces());
+
+ std::vector<Point> 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];
+
+ 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))
+ {
+ 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);
+ if (useAgglomeration)
+ {
+ startAgg.append
+ (
+ globalNumbering.toGlobal(proci, fAgg)
+ );
+ endAgg.append
+ (
+ globalNumbering.toGlobal(proci, remAgg)
+ );
+ }
+
+ label globalI = globalNumbering.toGlobal(proci, j);
+ endIndex.append(globalI);
+
+ 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 (debug)
+ {
+ Pout<< "Number of rays :" << totalRays << endl;
+ }
+
+ for (unsigned long rayI = 0; rayI < start.size(); ++rayI)
+ {
+ Segment ray(start[rayI], end[rayI]);
+ Segment_intersection intersects = tree.any_intersection(ray);
+
+ if (!intersects)
+ {
+ rayStartFace.append(startIndex[rayI]);
+ rayEndFace.append(endIndex[rayI]);
+ }
+ }
+ start.clear();
+
+ if (debug)
+ {
+ Pout << "hits : " << rayStartFace.size()<< endl;
+ }
+
+ startIndex.clear();
+ end.clear();
+ endIndex.clear();
+ startAgg.clear();
+ endAgg.clear();
+}
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index 264f141cd20..d34a5eb0acb 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,7 +25,7 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
- viewFactorsGenExt
+ viewFactorsGen
Group
grpPreProcessingUtilities
@@ -33,7 +33,7 @@ Group
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'.
+ when the two surfaces are 'far' apart and 2LI when 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.
@@ -57,8 +57,8 @@ Description
dumpRays dumps rays
- In order to specify the participants patches in the VF calculation the
- keywaord viewFactorWall should be added to the boundary file.
+ The participating patches in the VF calculation have to be in the
+ 'viewFactorWall' patch group (in the polyMesh/boundary file), e.g.
floor
{
@@ -68,6 +68,9 @@ Description
startFace 3100;
}
+ Compile with -DNO_CGAL only if no CGAL present - CGAL AABB tree performs
+ better than the built-in octree.
+
\*---------------------------------------------------------------------------*/
#include "argList.H"
@@ -181,15 +184,13 @@ triSurface triangulate
//triSurfaceToAgglom.resize(localTriFaceI-1);
triangles.shrink();
+ triSurface surface(triangles, mesh.points());
+ surface.compactPoints();
- triSurface rawSurface(triangles, mesh.points());
- triSurface surface
- (
- rawSurface.localFaces(),
- rawSurface.localPoints()
- );
+#ifndef NO_CGAL
+ // CGAL : every processor has whole surface
globalIndex globalFaceIdx(surface.size(), globalIndex::gatherOnly());
globalIndex globalPointIdx
@@ -228,6 +229,7 @@ triSurface triangulate
);
Pstream::broadcast(surface);
+#endif
// Add patch names to surface
surface.patches().setSize(newPatchI);
@@ -337,6 +339,7 @@ void insertMatrixElements
}
}
+
scalar GaussQuad
(
@@ -403,6 +406,8 @@ scalar GaussQuad
}
return dIntFij;
}
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
@@ -643,7 +648,13 @@ int main(int argc, char *argv[])
// Set up searching engine for obstacles
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#ifdef NO_CGAL
+ // Using octree
#include "searchingEngine.H"
+#else
+ // Using CGAL aabbtree (faster, more robust)
+ #include "searchingEngine_CGAL.H"
+#endif
// Determine rays between coarse face centres
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -653,7 +664,13 @@ int main(int argc, char *argv[])
// Return rayStartFace in local index and rayEndFace in global index
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#ifdef NO_CGAL
+ // Using octree, distributedTriSurfaceMesh
#include "shootRays.H"
+#else
+ // Using CGAL aabbtree (faster, more robust)
+ #include "shootRays_CGAL.H"
+#endif
// Calculate number of visible faces from local index
labelList nVisibleFaceFaces(nCoarseFaces, Zero);
@@ -1135,7 +1152,7 @@ int main(int argc, char *argv[])
forAll(coarseToFine, coarseI)
{
- scalar FiSum = sum(F2LI[compactI]);
+ const scalar FiSum = sum(F2LI[compactI]);
const label coarseFaceID = coarsePatchFace[coarseI];
const labelList& fineFaces = coarseToFine[coarseFaceID];
--
GitLab