ENH: AMI - refactored partialFaceAreaWeightAMI to derive from faceAreaWeightAMI

This commit is contained in:
andy 2013-05-28 15:06:33 +01:00
parent 16b2973fab
commit 0aa5613dcf
4 changed files with 167 additions and 538 deletions

View File

@ -25,7 +25,85 @@ License
#include "faceAreaWeightAMI.H" #include "faceAreaWeightAMI.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
(
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
List<DynamicList<label> >& tgtAddr,
List<DynamicList<scalar> >& tgtWght,
label srcFaceI,
label tgtFaceI
)
{
// construct weights and addressing
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label nFacesRemaining = srcAddr.size();
// 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);
// list to keep track of tgt faces used to seed src faces
labelList seedFaces(nFacesRemaining, -1);
seedFaces[srcFaceI] = tgtFaceI;
// list to keep track of whether src face can be mapped
boolList mapFlag(nFacesRemaining, true);
// reset starting seed
label startSeedI = 0;
DynamicList<label> nonOverlapFaces;
do
{
// Do advancing front starting from srcFaceI,tgtFaceI
bool faceProcessed = processSourceFace
(
srcFaceI,
tgtFaceI,
nbrFaces,
visitedFaces,
srcAddr,
srcWght,
tgtAddr,
tgtWght
);
mapFlag[srcFaceI] = false;
nFacesRemaining--;
if (!faceProcessed)
{
nonOverlapFaces.append(srcFaceI);
}
// choose new src face from current src face neighbour
if (nFacesRemaining > 0)
{
setNextFaces
(
startSeedI,
srcFaceI,
tgtFaceI,
mapFlag,
seedFaces,
visitedFaces
);
}
} while (nFacesRemaining > 0);
this->srcNonOverlap_.transfer(nonOverlapFaces);
}
template<class SourcePatch, class TargetPatch> template<class SourcePatch, class TargetPatch>
bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
@ -45,6 +123,11 @@ bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
List<DynamicList<scalar> >& tgtWght List<DynamicList<scalar> >& tgtWght
) )
{ {
if (tgtStartFaceI == -1)
{
return false;
}
nbrFaces.clear(); nbrFaces.clear();
visitedFaces.clear(); visitedFaces.clear();
@ -101,11 +184,15 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
label& tgtFaceI, label& tgtFaceI,
const boolList& mapFlag, const boolList& mapFlag,
labelList& seedFaces, labelList& seedFaces,
const DynamicList<label>& visitedFaces const DynamicList<label>& visitedFaces,
bool errorOnNotFound
) const ) const
{ {
const labelList& srcNbrFaces = this->srcPatch_.faceFaces()[srcFaceI]; const labelList& srcNbrFaces = this->srcPatch_.faceFaces()[srcFaceI];
// initialise tgtFaceI
tgtFaceI = -1;
// set possible seeds for later use // set possible seeds for later use
bool valuesSet = false; bool valuesSet = false;
forAll(srcNbrFaces, i) forAll(srcNbrFaces, i)
@ -195,19 +282,23 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
} }
} }
FatalErrorIn if (errorOnNotFound)
( {
"void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::" FatalErrorIn
"setNextFaces" (
"(" "void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::"
"label&, " "setNextFaces"
"label&, " "("
"label&, " "label&, "
"const boolList&, " "label&, "
"labelList&, " "label&, "
"const DynamicList<label>&" "const boolList&, "
") const" "labelList&, "
) << "Unable to set source and target faces" << abort(FatalError); "const DynamicList<label>&, "
"bool"
") const"
) << "Unable to set source and target faces" << abort(FatalError);
}
} }
} }
@ -447,77 +538,21 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
List<DynamicList<label> > tgtAddr(this->tgtPatch_.size()); List<DynamicList<label> > tgtAddr(this->tgtPatch_.size());
List<DynamicList<scalar> > tgtWght(tgtAddr.size()); List<DynamicList<scalar> > tgtWght(tgtAddr.size());
calcAddressing
(
srcAddr,
srcWght,
tgtAddr,
tgtWght,
srcFaceI,
tgtFaceI
);
// construct weights and addressing if (this->srcNonOverlap_.size() != 0)
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label nFacesRemaining = srcAddr.size();
// 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);
// list to keep track of tgt faces used to seed src faces
labelList seedFaces(nFacesRemaining, -1);
seedFaces[srcFaceI] = tgtFaceI;
// list to keep track of whether src face can be mapped
boolList mapFlag(nFacesRemaining, true);
// reset starting seed
label startSeedI = 0;
DynamicList<label> nonOverlapFaces;
do
{ {
// Do advancing front starting from srcFaceI,tgtFaceI Pout<< " AMI: " << this->srcNonOverlap_.size()
bool faceProcessed = processSourceFace
(
srcFaceI,
tgtFaceI,
nbrFaces,
visitedFaces,
srcAddr,
srcWght,
tgtAddr,
tgtWght
);
mapFlag[srcFaceI] = false;
nFacesRemaining--;
if (!faceProcessed)
{
nonOverlapFaces.append(srcFaceI);
}
// choose new src face from current src face neighbour
if (nFacesRemaining > 0)
{
setNextFaces
(
startSeedI,
srcFaceI,
tgtFaceI,
mapFlag,
seedFaces,
visitedFaces
);
}
} while (nFacesRemaining > 0);
if (nonOverlapFaces.size() != 0)
{
Pout<< " AMI: " << nonOverlapFaces.size()
<< " non-overlap faces identified" << " non-overlap faces identified"
<< endl; << endl;
this->srcNonOverlap_.transfer(nonOverlapFaces);
} }

View File

@ -52,9 +52,9 @@ class faceAreaWeightAMI
public AMIMethod<SourcePatch, TargetPatch> public AMIMethod<SourcePatch, TargetPatch>
{ {
private: protected:
// Private Member Functions // Protected Member Functions
//- Disallow default bitwise copy construct //- Disallow default bitwise copy construct
faceAreaWeightAMI(const faceAreaWeightAMI&); faceAreaWeightAMI(const faceAreaWeightAMI&);
@ -64,8 +64,19 @@ private:
// Marching front // Marching front
//- Calculate addressing and weights using temporary storage
virtual void calcAddressing
(
List<DynamicList<label> >& srcAddress,
List<DynamicList<scalar> >& srcWeights,
List<DynamicList<label> >& tgtAddress,
List<DynamicList<scalar> >& tgtWeights,
label srcFaceI,
label tgtFaceI
);
//- Determine overlap contributions for source face srcFaceI //- Determine overlap contributions for source face srcFaceI
bool processSourceFace virtual bool processSourceFace
( (
const label srcFaceI, const label srcFaceI,
const label tgtStartFaceI, const label tgtStartFaceI,
@ -78,7 +89,7 @@ private:
); );
//- Attempt to re-evaluate source faces that have not been included //- Attempt to re-evaluate source faces that have not been included
void restartUncoveredSourceFace virtual void restartUncoveredSourceFace
( (
List<DynamicList<label> >& srcAddr, List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght, List<DynamicList<scalar> >& srcWght,
@ -87,21 +98,22 @@ private:
); );
//- Set the source and target seed faces //- Set the source and target seed faces
void setNextFaces virtual void setNextFaces
( (
label& startSeedI, label& startSeedI,
label& srcFaceI, label& srcFaceI,
label& tgtFaceI, label& tgtFaceI,
const boolList& mapFlag, const boolList& mapFlag,
labelList& seedFaces, labelList& seedFaces,
const DynamicList<label>& visitedFaces const DynamicList<label>& visitedFaces,
bool errorOnNotFound = true
) const; ) const;
// Evaluation // Evaluation
//- Area of intersection between source and target faces //- Area of intersection between source and target faces
scalar interArea virtual scalar interArea
( (
const label srcFaceI, const label srcFaceI,
const label tgtFaceI const label tgtFaceI

View File

@ -27,77 +27,6 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
bool Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
(
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
)
{
if (tgtStartFaceI == -1)
{
return false;
}
nbrFaces.clear();
visitedFaces.clear();
// append initial target face and neighbours
nbrFaces.append(tgtStartFaceI);
this->appendNbrFaces
(
tgtStartFaceI,
this->tgtPatch_,
visitedFaces,
nbrFaces
);
bool faceProcessed = false;
do
{
// process new target face
label tgtFaceI = nbrFaces.remove();
visitedFaces.append(tgtFaceI);
scalar area = interArea(srcFaceI, tgtFaceI);
// store when intersection fractional area > tolerance
if (area/this->srcMagSf_[srcFaceI] > faceAreaIntersect::tolerance())
{
srcAddr[srcFaceI].append(tgtFaceI);
srcWght[srcFaceI].append(area);
tgtAddr[tgtFaceI].append(srcFaceI);
tgtWght[tgtFaceI].append(area);
this->appendNbrFaces
(
tgtFaceI,
this->tgtPatch_,
visitedFaces,
nbrFaces
);
faceProcessed = true;
}
} while (nbrFaces.size() > 0);
return faceProcessed;
}
template<class SourcePatch, class TargetPatch> template<class SourcePatch, class TargetPatch>
void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
( (
@ -106,272 +35,20 @@ void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
label& tgtFaceI, label& tgtFaceI,
const boolList& mapFlag, const boolList& mapFlag,
labelList& seedFaces, labelList& seedFaces,
const DynamicList<label>& visitedFaces const DynamicList<label>& visitedFaces,
const bool errorOnNotFound
) const ) const
{ {
const labelList& srcNbrFaces = this->srcPatch_.faceFaces()[srcFaceI]; faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
(
// set possible seeds for later use startSeedI,
bool valuesSet = false; srcFaceI,
forAll(srcNbrFaces, i) tgtFaceI,
{ mapFlag,
label faceS = srcNbrFaces[i]; seedFaces,
visitedFaces,
if (mapFlag[faceS] && seedFaces[faceS] == -1) false // no error on not found
{ );
forAll(visitedFaces, j)
{
label faceT = visitedFaces[j];
scalar area = interArea(faceS, faceT);
scalar areaTotal = this->srcMagSf_[srcFaceI];
// Check that faces have enough overlap for robust walking
if (area/areaTotal > faceAreaIntersect::tolerance())
{
// TODO - throwing area away - re-use in next iteration?
seedFaces[faceS] = faceT;
if (!valuesSet)
{
srcFaceI = faceS;
tgtFaceI = faceT;
valuesSet = true;
}
}
}
}
}
// set next src and tgt faces if not set above
if (valuesSet)
{
return;
}
else
{
// try to use existing seed
bool foundNextSeed = false;
for (label faceI = startSeedI; faceI < mapFlag.size(); faceI++)
{
if (mapFlag[faceI])
{
if (!foundNextSeed)
{
startSeedI = faceI;
foundNextSeed = true;
}
if (seedFaces[faceI] != -1)
{
srcFaceI = faceI;
tgtFaceI = seedFaces[faceI];
return;
}
}
}
// perform new search to find match
if (debug)
{
Pout<< "Advancing front stalled: searching for new "
<< "target face" << endl;
}
foundNextSeed = false;
for (label faceI = startSeedI; faceI < mapFlag.size(); faceI++)
{
if (mapFlag[faceI])
{
if (!foundNextSeed)
{
startSeedI = faceI + 1;
foundNextSeed = true;
}
srcFaceI = faceI;
tgtFaceI = this->findTargetFace(srcFaceI);
if (tgtFaceI >= 0)
{
return;
}
}
}
// no more matches to be found
tgtFaceI = -1;
return;
}
}
template<class SourcePatch, class TargetPatch>
Foam::scalar Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
(
const label srcFaceI,
const label tgtFaceI
) const
{
const pointField& srcPoints = this->srcPatch_.points();
const pointField& tgtPoints = this->tgtPatch_.points();
// references to candidate faces
const face& src = this->srcPatch_[srcFaceI];
const face& tgt = this->tgtPatch_[tgtFaceI];
// quick reject if either face has zero area
// Note: do not use stored face areas for target patch
const scalar tgtMag = tgt.mag(tgtPoints);
if ((this->srcMagSf_[srcFaceI] < ROOTVSMALL) || (tgtMag < ROOTVSMALL))
{
return 0.0;
}
// create intersection object
faceAreaIntersect inter(srcPoints, tgtPoints, this->reverseTarget_);
// crude resultant norm
vector n(-src.normal(srcPoints));
n /= mag(n);
if (this->reverseTarget_)
{
n -= tgt.normal(tgtPoints)/tgtMag;
}
else
{
n += tgt.normal(tgtPoints)/tgtMag;
}
n *= 0.5;
scalar area = 0;
if (mag(n) > ROOTVSMALL)
{
area = inter.calc(src, tgt, n, this->triMode_);
}
else
{
WarningIn
(
"void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::"
"interArea"
"("
"const label, "
"const label"
") 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))
{
this->writeIntersectionOBJ(area, src, tgt, srcPoints, tgtPoints);
}
return area;
}
template<class SourcePatch, class TargetPatch>
void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::
restartUncoveredSourceFace
(
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/this->srcMagSf_[srcFaceI];
if (t < 0.5)
{
lowWeightFaces.insert(srcFaceI);
}
}
if (debug)
{
Pout<< "partialFaceAreaWeightAMI: restarting search on "
<< lowWeightFaces.size() << " 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 = this->findTargetFace(srcFaceI);
if (tgtFaceI != -1)
{
//bool faceProcessed =
processSourceFace
(
srcFaceI,
tgtFaceI,
nbrFaces,
visitedFaces,
srcAddr,
srcWght,
tgtAddr,
tgtWght
);
// ? Check faceProcessed to see if restarting has worked.
}
}
}
} }
@ -389,7 +66,7 @@ partialFaceAreaWeightAMI
const bool reverseTarget const bool reverseTarget
) )
: :
AMIMethod<SourcePatch, TargetPatch> faceAreaWeightAMI<SourcePatch, TargetPatch>
( (
srcPatch, srcPatch,
tgtPatch, tgtPatch,
@ -412,7 +89,7 @@ Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch> template<class SourcePatch, class TargetPatch>
bool Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::normalise() const bool Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::conformal() const
{ {
return false; return false;
} }
@ -451,79 +128,15 @@ void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
List<DynamicList<label> > tgtAddr(this->tgtPatch_.size()); List<DynamicList<label> > tgtAddr(this->tgtPatch_.size());
List<DynamicList<scalar> > tgtWght(tgtAddr.size()); List<DynamicList<scalar> > tgtWght(tgtAddr.size());
faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
// construct weights and addressing (
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ srcAddr,
srcWght,
label nFacesRemaining = srcAddr.size(); tgtAddr,
tgtWght,
// list of tgt face neighbour faces srcFaceI,
DynamicList<label> nbrFaces(10); tgtFaceI
);
// list of faces currently visited for srcFaceI to avoid multiple hits
DynamicList<label> visitedFaces(10);
// list to keep track of tgt faces used to seed src faces
labelList seedFaces(nFacesRemaining, -1);
seedFaces[srcFaceI] = tgtFaceI;
// list to keep track of whether src face can be mapped
boolList mapFlag(nFacesRemaining, true);
// reset starting seed
label startSeedI = 0;
DynamicList<label> nonOverlapFaces;
do
{
// Do advancing front starting from srcFaceI,tgtFaceI
bool faceProcessed = processSourceFace
(
srcFaceI,
tgtFaceI,
nbrFaces,
visitedFaces,
srcAddr,
srcWght,
tgtAddr,
tgtWght
);
mapFlag[srcFaceI] = false;
nFacesRemaining--;
if (!faceProcessed)
{
nonOverlapFaces.append(srcFaceI);
}
// choose new src face from current src face neighbour
if (nFacesRemaining > 0)
{
setNextFaces
(
startSeedI,
srcFaceI,
tgtFaceI,
mapFlag,
seedFaces,
visitedFaces
);
}
} while (nFacesRemaining > 0);
if (nonOverlapFaces.size() != 0)
{
Pout<< " AMI: " << nonOverlapFaces.size()
<< " non-overlap faces identified"
<< endl;
this->srcNonOverlap_.transfer(nonOverlapFaces);
}
// transfer data to persistent storage // transfer data to persistent storage
forAll(srcAddr, i) forAll(srcAddr, i)

View File

@ -25,7 +25,7 @@ Class
Foam::partialFaceAreaWeightAMI Foam::partialFaceAreaWeightAMI
Description Description
Face area weighted Arbitrary Mesh Interface (AMI) method Partial face area weighted Arbitrary Mesh Interface (AMI) method
SourceFiles SourceFiles
partialFaceAreaWeightAMI.C partialFaceAreaWeightAMI.C
@ -35,7 +35,7 @@ SourceFiles
#ifndef partialFaceAreaWeightAMI_H #ifndef partialFaceAreaWeightAMI_H
#define partialFaceAreaWeightAMI_H #define partialFaceAreaWeightAMI_H
#include "AMIMethod.H" #include "faceAreaWeightAMI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -49,7 +49,7 @@ namespace Foam
template<class SourcePatch, class TargetPatch> template<class SourcePatch, class TargetPatch>
class partialFaceAreaWeightAMI class partialFaceAreaWeightAMI
: :
public AMIMethod<SourcePatch, TargetPatch> public faceAreaWeightAMI<SourcePatch, TargetPatch>
{ {
private: private:
@ -64,47 +64,16 @@ private:
// Marching front // Marching front
//- Determine overlap contributions for source face srcFaceI
bool processSourceFace
(
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
);
//- Attempt to re-evaluate source faces that have not been included
void restartUncoveredSourceFace
(
List<DynamicList<label> >& srcAddr,
List<DynamicList<scalar> >& srcWght,
List<DynamicList<label> >& tgtAddr,
List<DynamicList<scalar> >& tgtWght
);
//- Set the source and target seed faces //- Set the source and target seed faces
void setNextFaces virtual void setNextFaces
( (
label& startSeedI, label& startSeedI,
label& srcFaceI, label& srcFaceI,
label& tgtFaceI, label& tgtFaceI,
const boolList& mapFlag, const boolList& mapFlag,
labelList& seedFaces, labelList& seedFaces,
const DynamicList<label>& visitedFaces const DynamicList<label>& visitedFaces,
) const; bool errorOnNotFound = true
// Evaluation
//- Area of intersection between source and target faces
scalar interArea
(
const label srcFaceI,
const label tgtFaceI
) const; ) const;
@ -136,8 +105,8 @@ public:
// Access // Access
//- Flag to indicate that weights should be normalised //- Flag to indicate that interpolation patches are conformal
virtual bool normalise() const; virtual bool conformal() const;
// Manipulation // Manipulation