Merge branch 'feature-patchSeedSet' into 'develop'

Extended patchSeedSet to generate seed points more uniformly

See merge request Development/openfoam!647
This commit is contained in:
Kutalmış Berçin
2023-12-08 11:42:22 +00:00
5 changed files with 527 additions and 39 deletions

View File

@ -352,4 +352,6 @@ regionModel/regionProperties/regionProperties.C
tetOverlapVolume/tetOverlapVolume.C
triangulatedPatch/triangulatedPatch.C
LIB = $(FOAM_LIBBIN)/libmeshTools

View File

@ -0,0 +1,224 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023 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/>.
\*---------------------------------------------------------------------------*/
#include "triangulatedPatch.H"
#include "triPointRef.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::triangulatedPatch::randomPoint
(
Random& rnd,
const scalar c,
point& result,
label& facei,
label& celli
) const
{
result = point::min;
facei = -1;
celli = -1;
if (triWght_.empty() || c < triWght_.front() || c > triWght_.back())
{
return false;
}
// Find corresponding decomposed face triangle
// Note: triWght_ is sized nTri+1 (zero added at start)
label trii = 0;
for (label i = 0; i < triWght_.size() - 1; ++i)
{
if (c > triWght_[i] && c <= triWght_[i+1])
{
trii = i;
break;
}
}
// Find random point in triangle
const pointField& points = patch_.points();
const face& tf = triFace_[trii];
const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
result = tri.randomPoint(rnd);
facei = triToFace_[trii];
celli = patch_.faceCells()[facei];
if (perturbTol_ > 0)
{
const polyMesh& mesh = patch_.boundaryMesh().mesh();
const point& cc = mesh.cellCentres()[celli];
// Normal points out of domain => subtract correction
const vector& n = patch_.faceNormals()[facei];
result -= perturbTol_*n*mag(n & (cc - result));
// Reset facei - point no longer resides on the face
facei = -1;
}
return true;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::triangulatedPatch::triangulatedPatch
(
const polyPatch& patch,
const scalar perturbTol
)
:
patch_(patch),
perturbTol_(perturbTol),
triFace_(),
triToFace_(),
triWght_()
{
update();
}
Foam::triangulatedPatch::triangulatedPatch
(
const polyMesh& mesh,
const word& patchName,
const scalar perturbTol
)
:
triangulatedPatch(mesh.boundaryMesh()[patchName], perturbTol)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::triangulatedPatch::update()
{
const pointField& points = patch_.points();
// Triangulate the patch faces and create addressing
DynamicList<label> triToFace(2*patch_.size());
DynamicList<face> triFace(2*patch_.size());
DynamicList<scalar> triWght(2*patch_.size());
DynamicList<face> tris(8);
// Set zero value at the start of the tri area/weight list
triWght.push_back(0);
forAll(patch_, facei)
{
const face& f = patch_[facei];
tris.clear();
f.triangles(points, tris);
for (const auto& t : tris)
{
triToFace.push_back(facei);
triFace.push_back(t);
triWght.push_back(t.mag(points));
}
}
scalarList procSumWght(Pstream::nProcs()+1, Zero);
procSumWght[Pstream::myProcNo()+1] = sum(triWght);
Pstream::listCombineReduce(procSumWght, maxEqOp<scalar>());
for (label i = 1; i < procSumWght.size(); ++i)
{
// Convert to cumulative
procSumWght[i] += procSumWght[i-1];
}
const scalar offset = procSumWght[Pstream::myProcNo()];
forAll(triWght, i)
{
if (i)
{
// Convert to cumulative
triWght[i] += triWght[i-1];
}
// Apply processor offset
triWght[i] += offset;
}
// Normalise
const scalar sumWght = procSumWght.back();
for (scalar& w : triWght)
{
w /= sumWght;
}
// Transfer to persistent storage
triFace_.transfer(triFace);
triToFace_.transfer(triToFace);
triWght_.transfer(triWght);
}
bool Foam::triangulatedPatch::randomLocalPoint
(
Random& rnd,
point& result,
label& facei,
label& celli
) const
{
const scalar c = rnd.position<scalar>(triWght_.front(), triWght_.back());
return randomPoint(rnd, c, result, facei, celli);
}
bool Foam::triangulatedPatch::randomGlobalPoint
(
Random& rnd,
point& result,
label& facei,
label& celli
) const
{
boolList valid(UPstream::nProcs(), false);
valid[UPstream::myProcNo()] = randomLocalPoint(rnd, result, facei, celli);
UPstream::listGatherValues(valid);
forAll(valid, proci)
{
// Choose first valid processor
if (valid[proci])
{
return (proci == UPstream::myProcNo());
}
}
return false;
}
// ************************************************************************* //

View File

@ -0,0 +1,152 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023 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/>.
Class
Foam::triangulatedPatch
Description
Performs a triangulation of a patch to return randomised point locations.
\*---------------------------------------------------------------------------*/
#ifndef Foam_triangulatedPatch_H
#define Foam_triangulatedPatch_H
#include "polyMesh.H"
#include "Random.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class triangulatedPatch Declaration
\*---------------------------------------------------------------------------*/
class triangulatedPatch
{
// Private Data
//- Reference to the polyPatch
const polyPatch& patch_;
//- Perturbation tolerance to move the point towards the cell centre
bool perturbTol_;
//- Face triangles
faceList triFace_;
//- Triangle to patch face addressing
labelList triToFace_;
//- Triangle weights
scalarList triWght_;
// Private Member Functions
//- Update triangulation
void update();
//- Return a random point on the patch based on random number, 0 < c < 1
bool randomPoint
(
Random& rnd,
const scalar c,
point& result,
label& facei,
label& celli
) const;
public:
//- Constructors
//- Construct from components
triangulatedPatch
(
const polyPatch& patch,
const scalar perturbTol
);
//- Construct from mesh and patch name
triangulatedPatch
(
const polyMesh& mesh,
const word& patchName,
const scalar perturbTol
);
//- Destructor
~triangulatedPatch() = default;
// Member Functions
//- Set a random point on the local patch
//
// \param rnd random number generator
// \param result the random point
// \param facei index associated with the random point
// \param celli index associated with the random point
//
// \return true if point has been set
bool randomLocalPoint
(
Random& rnd,
point& result,
label& facei,
label& celli
) const;
//- Set a global random point on the patch
//
// \param rnd random number generator
// \param result the random point
// \param facei index associated with the random point
// \param celli index associated with the random point
//
// \return true if point has been set
bool randomGlobalPoint
(
Random& rnd,
point& result,
label& facei,
label& celli
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2012-2016 OpenFOAM Foundation
Copyright (C) 2018-2022 OpenCFD Ltd.
Copyright (C) 2018-2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -28,13 +28,12 @@ License
#include "patchSeedSet.H"
#include "polyMesh.H"
#include "addToRunTimeSelectionTable.H"
#include "treeBoundBox.H"
#include "treeDataFace.H"
#include "Time.H"
#include "meshTools.H"
#include "mappedPatchBase.H"
#include "indirectPrimitivePatch.H"
#include "triangulatedPatch.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -47,8 +46,10 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::patchSeedSet::calcSamples
void Foam::patchSeedSet::calcPatchSamples
(
const label nAvailable,
const label nPatchPoints,
DynamicList<point>& samplingPts,
DynamicList<label>& samplingCells,
DynamicList<label>& samplingFaces,
@ -56,39 +57,71 @@ void Foam::patchSeedSet::calcSamples
DynamicList<scalar>& samplingCurveDist
)
{
DebugInfo << "patchSeedSet : sampling on patches :" << endl;
if (nAvailable < 1)
{
return;
}
// Construct search tree for all patch faces.
Random& rndGen = *rndGenPtr_;
globalIndex globalSampleNumbers(nAvailable);
label nGlobalPatchPoints = returnReduce(nPatchPoints, sumOp<label>());
point pt;
label facei;
label celli;
const bool perturb = true;
for (const label patchi : patchSet_)
{
const polyPatch& pp = mesh().boundaryMesh()[patchi];
triangulatedPatch tp(pp, perturb);
const label np = nAvailable*pp.size()/scalar(nGlobalPatchPoints);
for (label i = 0; i < np; ++i)
{
tp.randomLocalPoint(rndGen, pt, facei, celli);
samplingPts.append(pt);
samplingCells.append(celli);
samplingFaces.append(pp.start() + facei);
samplingSegments.append(0);
samplingCurveDist.append(globalSampleNumbers.toGlobal(i));
}
}
}
void Foam::patchSeedSet::calcSelectedLocations
(
const label nAvailable,
const label nPatchPoints,
DynamicList<point>& samplingPts,
DynamicList<label>& samplingCells,
DynamicList<label>& samplingFaces,
DynamicList<label>& samplingSegments,
DynamicList<scalar>& samplingCurveDist
)
{
if (nAvailable < 1)
{
return;
}
Random& rndGen = *rndGenPtr_;
labelList patchFaces(nPatchPoints);
label sz = 0;
for (const label patchi : patchSet_)
{
const polyPatch& pp = mesh().boundaryMesh()[patchi];
sz += pp.size();
DebugInfo << " " << pp.name() << " size " << pp.size() << endl;
}
labelList patchFaces(sz);
sz = 0;
for (const label patchi : patchSet_)
{
const polyPatch& pp = mesh().boundaryMesh()[patchi];
forAll(pp, i)
forAll(pp, localFacei)
{
patchFaces[sz++] = pp.start()+i;
patchFaces[sz++] = pp.start() + localFacei;
}
}
if (!rndGenPtr_)
{
rndGenPtr_.reset(new Random(0));
}
Random& rndGen = *rndGenPtr_;
if (selectedLocations_.size())
{
DynamicList<label> newPatchFaces(patchFaces.size());
@ -98,7 +131,7 @@ void Foam::patchSeedSet::calcSamples
// selectedLocations
// All the info for nearest. Construct to miss
List<mappedPatchBase::nearInfo> nearest(selectedLocations_.size());
List<mappedPatchBase::nearInfo> nearest(nAvailable);
const indirectPrimitivePatch pp
(
@ -130,7 +163,7 @@ void Foam::patchSeedSet::calcSamples
GREAT
);
forAll(selectedLocations_, sampleI)
for (label sampleI = 0; sampleI < nAvailable; ++sampleI)
{
const auto& treeData = boundaryTree.shapes();
const point& sample = selectedLocations_[sampleI];
@ -184,7 +217,7 @@ void Foam::patchSeedSet::calcSamples
if (debug)
{
Pout<< "Found " << newPatchFaces.size()
<< " out of " << selectedLocations_.size()
<< " out of " << nAvailable
<< " on local processor" << endl;
}
@ -193,13 +226,12 @@ void Foam::patchSeedSet::calcSamples
// Shuffle and truncate if in random mode
label totalSize = returnReduce(patchFaces.size(), sumOp<label>());
const label totalSize = returnReduce(patchFaces.size(), sumOp<label>());
if (maxPoints_ < totalSize)
if (totalSize > nAvailable)
{
// Check what fraction of maxPoints_ I need to generate locally.
label myMaxPoints =
label(scalar(patchFaces.size())/totalSize*maxPoints_);
label myMaxPoints = scalar(patchFaces.size())/totalSize*nAvailable;
labelList subset = identity(patchFaces.size());
for (label iter = 0; iter < 4; ++iter)
@ -210,6 +242,7 @@ void Foam::patchSeedSet::calcSamples
std::swap(subset[i], subset[j]);
}
}
// Truncate
subset.setSize(myMaxPoints);
@ -239,7 +272,7 @@ void Foam::patchSeedSet::calcSamples
forAll(patchFaces, i)
{
label facei = patchFaces[i];
const label facei = patchFaces[i];
// Slightly shift point in since on warped face face-diagonal
// decomposition might be outside cell for face-centre decomposition!
@ -249,7 +282,7 @@ void Foam::patchSeedSet::calcSamples
facei,
polyMesh::FACE_DIAG_TRIS
);
label celli = mesh().faceOwner()[facei];
const label celli = mesh().faceOwner()[facei];
if (info.hit())
{
@ -313,6 +346,59 @@ void Foam::patchSeedSet::genSamples()
}
void Foam::patchSeedSet::calcSamples
(
DynamicList<point>& samplingPts,
DynamicList<label>& samplingCells,
DynamicList<label>& samplingFaces,
DynamicList<label>& samplingSegments,
DynamicList<scalar>& samplingCurveDist
)
{
DebugInfo << "patchSeedSet : sampling on patches :" << endl;
if (!rndGenPtr_)
{
rndGenPtr_.reset(new Random(0));
}
label nPatchPoints = 0;
for (const label patchi : patchSet_)
{
const polyPatch& pp = mesh().boundaryMesh()[patchi];
nPatchPoints += pp.size();
DebugInfo << " " << pp.name() << " size " << pp.size() << endl;
}
label nAvailable = min(maxPoints_, selectedLocations_.size());
calcSelectedLocations
(
nAvailable,
nPatchPoints,
samplingPts,
samplingCells,
samplingFaces,
samplingSegments,
samplingCurveDist
);
nAvailable = maxPoints_ - nAvailable;
calcPatchSamples
(
nAvailable,
nPatchPoints,
samplingPts,
samplingCells,
samplingFaces,
samplingSegments,
samplingCurveDist
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::patchSeedSet::patchSeedSet

View File

@ -58,7 +58,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class patchSeedSet Declaration
Class patchSeedSet Declaration
\*---------------------------------------------------------------------------*/
class patchSeedSet
@ -82,6 +82,30 @@ class patchSeedSet
// Private Member Functions
//- Determine (random) samples on patch
void calcPatchSamples
(
const label nAvailable,
const label nPatchPoints,
DynamicList<point>& samplingPts,
DynamicList<label>& samplingCells,
DynamicList<label>& samplingFaces,
DynamicList<label>& samplingSegments,
DynamicList<scalar>& samplingCurveDist
);
//- Determine samples provided by supplied points
void calcSelectedLocations
(
const label nAvailable,
const label nPatchPoints,
DynamicList<point>& samplingPts,
DynamicList<label>& samplingCells,
DynamicList<label>& samplingFaces,
DynamicList<label>& samplingSegments,
DynamicList<scalar>& samplingCurveDist
);
//- Samples all points in sampleCoords.
void calcSamples
(