ENH: AMIInterpolation: use stable guess for staritng seed. Add restart (under debug flag)

This commit is contained in:
mattijs
2012-07-25 12:09:09 +01:00
parent 49dc91fe4c
commit 83b357ab6e
5 changed files with 260 additions and 38 deletions

View File

@ -782,6 +782,62 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::appendNbrFaces
}
template<class SourcePatch, class TargetPatch>
bool Foam::AMIInterpolation<SourcePatch, TargetPatch>::processSourceFace
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const label srcFaceI,
const label tgtStartFaceI,
// list of tgt face neighbour faces
DynamicList<label>& nbrFaces,
// list of faces currently visited for srcFaceI to avoid multiple hits
DynamicList<label>& visitedFaces,
// temporary storage for addressing and weights
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
List<DynamicList<label> >& tgtAddr,
List<DynamicList<scalar> >& tgtWght
)
{
nbrFaces.clear();
visitedFaces.clear();
// append initial target face and neighbours
nbrFaces.append(tgtStartFaceI);
appendNbrFaces(tgtStartFaceI, tgtPatch, visitedFaces, nbrFaces);
bool faceProcessed = false;
do
{
// process new target face
label tgtFaceI = nbrFaces.remove();
visitedFaces.append(tgtFaceI);
scalar area = interArea(srcFaceI, tgtFaceI, srcPatch, tgtPatch);
// store when intersection area > 0
if (area > 0)
{
srcAddr[srcFaceI].append(tgtFaceI);
srcWght[srcFaceI].append(area);
tgtAddr[tgtFaceI].append(srcFaceI);
tgtWght[tgtFaceI].append(area);
appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
faceProcessed = true;
}
} while (nbrFaces.size() > 0);
return faceProcessed;
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
(
@ -811,7 +867,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
label faceT = visitedFaces[j];
scalar area = interArea(faceS, faceT, srcPatch0, tgtPatch0);
if (area > 0)
// Check that faces have enough overlap for robust walking
if (area/srcMagSf_[srcFaceI] > faceAreaIntersect::tolerance())
{
// TODO - throwing area away - re-use in next iteration?
@ -864,7 +921,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
<< "target face" << endl;
}
// foundNextSeed = false;
foundNextSeed = false;
for (label faceI = startSeedI; faceI < mapFlag.size(); faceI++)
{
if (mapFlag[faceI])
@ -887,7 +944,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
FatalErrorIn
(
"void Foam::cyclicAMIPolyPatch::setNextFaces"
"void Foam::AMIInterpolation<SourcePatch, TargetPatch>::"
"setNextFaces"
"("
"label&, "
"label&, "
@ -946,6 +1004,25 @@ Foam::scalar Foam::AMIInterpolation<SourcePatch, TargetPatch>::interArea
{
area = inter.calc(src, tgt, n, triMode_);
}
else
{
WarningIn
(
"void Foam::AMIInterpolation<SourcePatch, TargetPatch>::"
"interArea"
"("
"const label, "
"const label, "
"const SourcePatch&, "
"const TargetPatch&"
") const"
) << "Invalid normal for source face " << srcFaceI
<< " points " << UIndirectList<point>(srcPoints, src)
<< " target face " << tgtFaceI
<< " points " << UIndirectList<point>(tgtPoints, tgt)
<< endl;
}
if ((debug > 1) && (area > 0))
{
@ -956,6 +1033,105 @@ Foam::scalar Foam::AMIInterpolation<SourcePatch, TargetPatch>::interArea
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::
restartUncoveredSourceFace
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
List<DynamicList<label> >& tgtAddr,
List<DynamicList<scalar> >& tgtWght
)
{
// Collect all src faces with a low weight
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelHashSet lowWeightFaces(100);
forAll(srcWght, srcFaceI)
{
scalar s = sum(srcWght[srcFaceI]);
scalar t = s/srcMagSf_[srcFaceI];
if (t < 0.5)
{
lowWeightFaces.insert(srcFaceI);
}
}
Info<< "AMIInterpolation : restarting search on "
<< returnReduce(lowWeightFaces.size(), sumOp<label>())
<< " faces since sum of weights < 0.5" << endl;
if (lowWeightFaces.size() > 0)
{
// Erase all the lowWeight source faces from the target
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
DynamicList<label> okSrcFaces(10);
DynamicList<scalar> okSrcWeights(10);
forAll(tgtAddr, tgtFaceI)
{
okSrcFaces.clear();
okSrcWeights.clear();
DynamicList<label>& srcFaces = tgtAddr[tgtFaceI];
DynamicList<scalar>& srcWeights = tgtWght[tgtFaceI];
forAll(srcFaces, i)
{
if (!lowWeightFaces.found(srcFaces[i]))
{
okSrcFaces.append(srcFaces[i]);
okSrcWeights.append(srcWeights[i]);
}
}
if (okSrcFaces.size() < srcFaces.size())
{
srcFaces.transfer(okSrcFaces);
srcWeights.transfer(okSrcWeights);
}
}
// Restart search from best hit
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// list of tgt face neighbour faces
DynamicList<label> nbrFaces(10);
// list of faces currently visited for srcFaceI to avoid multiple hits
DynamicList<label> visitedFaces(10);
forAllConstIter(labelHashSet, lowWeightFaces, iter)
{
label srcFaceI = iter.key();
label tgtFaceI = findTargetFace(srcFaceI, srcPatch);
if (tgtFaceI != -1)
{
//bool faceProcessed =
processSourceFace
(
srcPatch,
tgtPatch,
srcFaceI,
tgtFaceI,
nbrFaces,
visitedFaces,
srcAddr,
srcWght,
tgtAddr,
tgtWght
);
// ? Check faceProcessed to see if restarting has worked.
}
}
}
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
(
@ -1073,37 +1249,22 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
DynamicList<label> nonOverlapFaces;
do
{
nbrFaces.clear();
visitedFaces.clear();
// Do advancing front starting from srcFaceI,tgtFaceI
bool faceProcessed = processSourceFace
(
srcPatch,
tgtPatch,
srcFaceI,
tgtFaceI,
// append initial target face and neighbours
nbrFaces.append(tgtFaceI);
appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
nbrFaces,
visitedFaces,
bool faceProcessed = false;
do
{
// process new target face
tgtFaceI = nbrFaces.remove();
visitedFaces.append(tgtFaceI);
scalar area = interArea(srcFaceI, tgtFaceI, srcPatch, tgtPatch);
// store when intersection area > 0
if (area > 0)
{
srcAddr[srcFaceI].append(tgtFaceI);
srcWght[srcFaceI].append(area);
tgtAddr[tgtFaceI].append(srcFaceI);
tgtWght[tgtFaceI].append(area);
appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
faceProcessed = true;
}
} while (nbrFaces.size() > 0);
srcAddr,
srcWght,
tgtAddr,
tgtWght
);
mapFlag[srcFaceI] = false;
@ -1140,6 +1301,22 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
srcNonOverlap_.transfer(nonOverlapFaces);
}
// Check for badly covered faces
if (debug)
{
restartUncoveredSourceFace
(
srcPatch,
tgtPatch,
srcAddr,
srcWght,
tgtAddr,
tgtWght
);
}
// transfer data to persistent storage
forAll(srcAddr, i)
{

View File

@ -246,6 +246,31 @@ class AMIInterpolation
DynamicList<label>& faceIDs
) const;
bool processSourceFace
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const label srcFaceI,
const label tgtStartFaceI,
DynamicList<label>& nbrFaces,
DynamicList<label>& visitedFaces,
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
List<DynamicList<label> >& tgtAddr,
List<DynamicList<scalar> >& tgtWght
);
void restartUncoveredSourceFace
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
List<DynamicList<label> >& tgtAddr,
List<DynamicList<scalar> >& tgtWght
);
//- Set the source and target seed faces
void setNextFaces
(

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,6 +25,10 @@ License
#include "faceAreaIntersect.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
Foam::scalar Foam::faceAreaIntersect::tol = 1e-6;
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
void Foam::faceAreaIntersect::triSliceWithPlane
@ -36,8 +40,6 @@ void Foam::faceAreaIntersect::triSliceWithPlane
const scalar len
)
{
const scalar matchTol = 1e-6;
// distance to cutting plane
FixedList<scalar, 3> d;
@ -51,7 +53,7 @@ void Foam::faceAreaIntersect::triSliceWithPlane
{
d[i] = ((tri[i] - p.refPoint()) & p.normal());
if (mag(d[i]) < matchTol*len)
if (mag(d[i]) < tol*len)
{
nCoPlanar++;
copI = i;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -78,6 +78,11 @@ private:
const bool reverseB_;
// Static data members
static scalar tol;
// Private Member Functions
//- Get triPoints from face
@ -152,6 +157,9 @@ public:
// Public Member Functions
//- Fraction of local length scale to use as intersection tolerance
inline static scalar& tolerance();
//- Return area of intersection of faceA with faceB
scalar calc
(

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,6 +23,8 @@ License
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
inline void Foam::faceAreaIntersect::setTriPoints
(
const point& a,
@ -106,4 +108,12 @@ inline Foam::scalar Foam::faceAreaIntersect::triArea(const triPoints& t) const
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::scalar& Foam::faceAreaIntersect::tolerance()
{
return tol;
}
// ************************************************************************* //