zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: viewFactorGen: revert to v2206 without CGAL

Browse Source
This commit is contained in:
mattijs
2022-08-04 14:24:50 +01:00
committed by Kutalmis Bercin
parent 994addc543
commit c08fc5ecd9
7 changed files with 371 additions and 85 deletions
Download Patch File Download Diff File Expand all files Collapse all files

19
applications/utilities/preProcessing/viewFactorsGen/Allwmake Executable file
View File

@ -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
#------------------------------------------------------------------------------

3
applications/utilities/preProcessing/viewFactorsGen/Make/options
View File

@ -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 \

31
applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
View File

@ -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());

84
applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H Normal file
View File

@ -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());

143
applications/utilities/preProcessing/viewFactorsGen/shootRays.H
View File

@ -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));
}
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)
// Return hit global agglo index
labelList aggHitIndex;
surfacesMesh.getField(hitInfo, aggHitIndex);
DynamicList<label> dRayIs;
// 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)
{
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)
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();

131
applications/utilities/preProcessing/viewFactorsGen/shootRays_CGAL.H Normal file
View File

@ -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();
}

41
applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
View File

@ -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];