ENH: AMI - added new partialFaceAreaWeightAMI option

This commit is contained in:
andy
2013-05-28 12:34:44 +01:00
parent 8205da653a
commit 2f8f797134
5 changed files with 790 additions and 50 deletions

View File

@ -56,6 +56,11 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolationMethodToWord
method = "faceAreaWeightAMI";
break;
}
case imPartialFaceAreaWeight:
{
method = "partialFaceAreaWeightAMI";
break;
}
default:
{
FatalErrorIn
@ -87,7 +92,15 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::wordTointerpolationMethod
wordList methods
(
IStringStream("(directAMI mapNearestAMI faceAreaWeightAMI)")()
IStringStream
(
"("
"directAMI "
"mapNearestAMI "
"faceAreaWeightAMI "
"partialFaceAreaWeightAMI"
")"
)()
);
if (im == "directAMI")
@ -102,6 +115,10 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::wordTointerpolationMethod
{
method = imFaceAreaWeight;
}
else if (im == "partialFaceAreaWeightAMI")
{
method = imPartialFaceAreaWeight;
}
else
{
FatalErrorIn
@ -183,6 +200,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights
const labelListList& addr,
scalarListList& wght,
scalarField& wghtSum,
const bool normalise,
const bool output
)
{
@ -191,12 +209,19 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights
forAll(wght, faceI)
{
scalarList& w = wght[faceI];
scalar denom = patchAreas[faceI];
scalar s = sum(w);
scalar t = s/patchAreas[faceI];
scalar t = s/denom;
if (normalise)
{
denom = s;
}
forAll(w, i)
{
w[i] /= s;
w[i] /= denom;
}
wghtSum[faceI] = t;
@ -476,6 +501,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
srcAddress,
srcWeights,
srcWeightsSum,
true,
false
);
}
@ -749,7 +775,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
tgtMagSf_[faceI] = tgtPatch[faceI].mag(tgtPatch.points());
}
// Calculate if patches present on multiple processors
singlePatchProc_ = calcDistribution(srcPatch, tgtPatch);
@ -794,9 +819,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
newTgtPoints
);
// calculate AMI interpolation
{
autoPtr<AMIMethod<SourcePatch, TargetPatch> > AMIPtr
(
AMIMethod<SourcePatch, TargetPatch>::New
@ -818,8 +841,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
tgtAddress_,
tgtWeights_
);
}
// Now
// ~~~
@ -886,6 +907,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
srcAddress_,
srcWeights_,
srcWeightsSum_,
AMIPtr->normalise(),
true
);
normaliseWeights
@ -895,6 +917,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
tgtAddress_,
tgtWeights_,
tgtWeightsSum_,
AMIPtr->normalise(),
true
);
@ -903,7 +926,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
srcMapPtr_.reset(new mapDistribute(globalSrcFaces, tgtAddress_, cMap));
tgtMapPtr_.reset(new mapDistribute(globalTgtFaces, srcAddress_, cMap));
if (debug)
{
writeFaceConnectivity(srcPatch, newTgtPatch, srcAddress_);
@ -913,7 +935,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
{
// calculate AMI interpolation
{
autoPtr<AMIMethod<SourcePatch, TargetPatch> > AMIPtr
(
AMIMethod<SourcePatch, TargetPatch>::New
@ -935,7 +956,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
tgtAddress_,
tgtWeights_
);
}
normaliseWeights
(
@ -944,6 +964,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
srcAddress_,
srcWeights_,
srcWeightsSum_,
AMIPtr->normalise(),
true
);
normaliseWeights
@ -953,6 +974,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
tgtAddress_,
tgtWeights_,
tgtWeightsSum_,
AMIPtr->normalise(),
true
);
}

View File

@ -88,7 +88,8 @@ public:
{
imDirect,
imMapNearest,
imFaceAreaWeight
imFaceAreaWeight,
imPartialFaceAreaWeight
};
//- Convert interpolationMethod to word representation
@ -236,6 +237,7 @@ private:
const labelListList& addr,
scalarListList& wght,
scalarField& wghtSum,
const bool normalise,
const bool output
);

View File

@ -0,0 +1,542 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 "partialFaceAreaWeightAMI.H"
// * * * * * * * * * * * * * 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>
void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
(
label& startSeedI,
label& srcFaceI,
label& tgtFaceI,
const boolList& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces
) const
{
const labelList& srcNbrFaces = this->srcPatch_.faceFaces()[srcFaceI];
// set possible seeds for later use
bool valuesSet = false;
forAll(srcNbrFaces, i)
{
label faceS = srcNbrFaces[i];
if (mapFlag[faceS] && seedFaces[faceS] == -1)
{
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.
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::
partialFaceAreaWeightAMI
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget
)
:
AMIMethod<SourcePatch, TargetPatch>
(
srcPatch,
tgtPatch,
srcMagSf,
tgtMagSf,
triMode,
reverseTarget
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::
~partialFaceAreaWeightAMI()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
bool Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::normalise() const
{
return false;
}
template<class SourcePatch, class TargetPatch>
void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFaceI,
label tgtFaceI
)
{
bool ok =
this->initialise
(
srcAddress,
srcWeights,
tgtAddress,
tgtWeights,
srcFaceI,
tgtFaceI
);
if (!ok)
{
return;
}
// temporary storage for addressing and weights
List<DynamicList<label> > srcAddr(this->srcPatch_.size());
List<DynamicList<scalar> > srcWght(srcAddr.size());
List<DynamicList<label> > tgtAddr(this->tgtPatch_.size());
List<DynamicList<scalar> > tgtWght(tgtAddr.size());
// 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);
if (nonOverlapFaces.size() != 0)
{
Pout<< " AMI: " << nonOverlapFaces.size()
<< " non-overlap faces identified"
<< endl;
this->srcNonOverlap_.transfer(nonOverlapFaces);
}
// transfer data to persistent storage
forAll(srcAddr, i)
{
srcAddress[i].transfer(srcAddr[i]);
srcWeights[i].transfer(srcWght[i]);
}
forAll(tgtAddr, i)
{
tgtAddress[i].transfer(tgtAddr[i]);
tgtWeights[i].transfer(tgtWght[i]);
}
}
// ************************************************************************* //

View File

@ -0,0 +1,172 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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::partialFaceAreaWeightAMI
Description
Face area weighted Arbitrary Mesh Interface (AMI) method
SourceFiles
partialFaceAreaWeightAMI.C
\*---------------------------------------------------------------------------*/
#ifndef partialFaceAreaWeightAMI_H
#define partialFaceAreaWeightAMI_H
#include "AMIMethod.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class partialFaceAreaWeightAMI Declaration
\*---------------------------------------------------------------------------*/
template<class SourcePatch, class TargetPatch>
class partialFaceAreaWeightAMI
:
public AMIMethod<SourcePatch, TargetPatch>
{
private:
// Private Member Functions
//- Disallow default bitwise copy construct
partialFaceAreaWeightAMI(const partialFaceAreaWeightAMI&);
//- Disallow default bitwise assignment
void operator=(const partialFaceAreaWeightAMI&);
// 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
void setNextFaces
(
label& startSeedI,
label& srcFaceI,
label& tgtFaceI,
const boolList& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces
) const;
// Evaluation
//- Area of intersection between source and target faces
scalar interArea
(
const label srcFaceI,
const label tgtFaceI
) const;
public:
//- Runtime type information
TypeName("partialFaceAreaWeightAMI");
// Constructors
//- Construct from components
partialFaceAreaWeightAMI
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch,
const scalarField& srcMagSf,
const scalarField& tgtMagSf,
const faceAreaIntersect::triangulationMode& triMode,
const bool reverseTarget = false
);
//- Destructor
virtual ~partialFaceAreaWeightAMI();
// Member Functions
// Access
//- Flag to indicate that weights should be normalised
virtual bool normalise() const;
// Manipulation
//- Update addressing and weights
virtual void calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFaceI = -1,
label tgtFaceI = -1
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "partialFaceAreaWeightAMI.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -28,6 +28,7 @@ License
#include "directAMI.H"
#include "mapNearestAMI.H"
#include "faceAreaWeightAMI.H"
#include "partialFaceAreaWeightAMI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -38,6 +39,7 @@ namespace Foam
makeAMIMethodType(AMIPatchToPatchInterpolation, directAMI);
makeAMIMethodType(AMIPatchToPatchInterpolation, mapNearestAMI);
makeAMIMethodType(AMIPatchToPatchInterpolation, faceAreaWeightAMI);
makeAMIMethodType(AMIPatchToPatchInterpolation, partialFaceAreaWeightAMI);
}