exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Updates to GGIInterpolation

Browse source 
Correct handling of partially overlapped faces. Initial commit, not tested.
This commit is contained in:
Vuko Vukcevic 2017-10-13 08:35:13 +02:00
parent f4db00045d
commit ba798d60ac
6 changed files with 301 additions and 18 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

24
src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C
View file

@ -211,18 +211,20 @@ void ggiFvPatchField<Type>::initEvaluate
+ (1.0 - this->patch().weights())*this->patchNeighbourField()
);
if (ggiPatch_.bridgeOverlap())
{
// Symmetry treatment used for overlap
vectorField nHat = this->patch().nf();
// Note: bridging already carried out when calling patchNeighbourField()
Field<Type> pif = this->patchInternalField();
Field<Type> bridgeField =
0.5*(pif + transform(I - 2.0*sqr(nHat), pif));
ggiPatch_.bridge(bridgeField, pf);
}
// if (ggiPatch_.bridgeOverlap())
// {
// // Symmetry treatment used for overlap
// vectorField nHat = this->patch().nf();
//
// Field<Type> pif = this->patchInternalField();
//
// Field<Type> bridgeField =
// 0.5*(pif + transform(I - 2.0*sqr(nHat), pif));
//
// ggiPatch_.bridge(bridgeField, pf);
// }
Field<Type>::operator=(pf);
}

31
src/foam/interpolations/GGIInterpolation/GGIInterpolate.C
View file

@ -97,13 +97,24 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridge
(
const Field<Type>& bridgeField,
Field<Type>& ff,
const labelList& addr
const labelList& addr,
const labelList& partiallyCoveredAddr,
const scalarField& uncoveredFractions
)
{
// Fully uncovered faces
forAll (addr, faceI)
{
ff[addr[faceI]] = bridgeField[addr[faceI]];
}
// Partially covered faces. Note the operator+= since we assume that the
// interpolation part is carried out before bridging (see
// e.g. ggiFvPatchField::patchNeighbourField())
forAll (partiallyCoveredAddr, faceI)
{
ff[addr[faceI]] += uncoveredFractions[faceI]*bridgeField[addr[faceI]];
}
}
@ -527,7 +538,14 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridgeMaster
<< abort(FatalError);
}
bridge(bridgeField, ff, uncoveredMasterFaces());
bridge
(
bridgeField,
ff,
uncoveredMasterFaces(),
partiallyCoveredMasterFaces(),
masterFaceUncoveredFractions()
);
}
@ -592,7 +610,14 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridgeSlave
<< abort(FatalError);
}
bridge(bridgeField, ff, uncoveredSlaveFaces());
bridge
(
bridgeField,
ff,
uncoveredSlaveFaces(),
partiallyCoveredSlaveFaces(),
slaveFaceUncoveredFractions()
);
}

62
src/foam/interpolations/GGIInterpolation/GGIInterpolation.C
View file

@ -104,7 +104,11 @@ void GGIInterpolation<MasterPatch, SlavePatch>::clearOut()
deleteDemandDrivenData(slaveWeightsPtr_);
deleteDemandDrivenData(uncoveredMasterAddrPtr_);
deleteDemandDrivenData(partiallyCoveredMasterAddrPtr_);
deleteDemandDrivenData(masterFaceUncoveredFractionsPtr_);
deleteDemandDrivenData(uncoveredSlaveAddrPtr_);
deleteDemandDrivenData(partiallyCoveredSlaveAddrPtr_);
deleteDemandDrivenData(slaveFaceUncoveredFractionsPtr_);
deleteDemandDrivenData(masterPointAddressingPtr_);
deleteDemandDrivenData(masterPointWeightsPtr_);
@ -155,7 +159,11 @@ GGIInterpolation<MasterPatch, SlavePatch>::GGIInterpolation
slavePointWeightsPtr_(NULL),
slavePointDistancePtr_(NULL),
uncoveredMasterAddrPtr_(NULL),
uncoveredSlaveAddrPtr_(NULL)
partiallyCoveredMasterAddrPtr_(NULL),
masterFaceUncoveredFractionsPtr_(NULL),
uncoveredSlaveAddrPtr_(NULL),
partiallyCoveredSlaveAddrPtr_(NULL),
slaveFaceUncoveredFractionsPtr_(NULL)
{
// Check size of transform. They should be equal to slave patch size
// if the transform is not constant
@ -256,6 +264,32 @@ GGIInterpolation<MasterPatch, SlavePatch>::uncoveredMasterFaces() const
}
template<class MasterPatch, class SlavePatch>
const labelList&
GGIInterpolation<MasterPatch, SlavePatch>::partiallyCoveredMasterFaces() const
{
if (!partiallyCoveredMasterAddrPtr_)
{
calcAddressing();
}
return *partiallyCoveredMasterAddrPtr_;
}
template<class MasterPatch, class SlavePatch>
const scalarField&
GGIInterpolation<MasterPatch, SlavePatch>::masterFaceUncoveredFractions() const
{
if (!masterFaceUncoveredFractionsPtr_)
{
calcAddressing();
}
return *masterFaceUncoveredFractionsPtr_;
}
template<class MasterPatch, class SlavePatch>
const labelList&
GGIInterpolation<MasterPatch, SlavePatch>::uncoveredSlaveFaces() const
@ -269,6 +303,32 @@ GGIInterpolation<MasterPatch, SlavePatch>::uncoveredSlaveFaces() const
}
template<class MasterPatch, class SlavePatch>
const labelList&
GGIInterpolation<MasterPatch, SlavePatch>::partiallyCoveredSlaveFaces() const
{
if (!partiallyCoveredSlaveAddrPtr_)
{
calcAddressing();
}
return *partiallyCoveredSlaveAddrPtr_;
}
template<class MasterPatch, class SlavePatch>
const scalarField&
GGIInterpolation<MasterPatch, SlavePatch>::slaveFaceUncoveredFractions() const
{
if (!slaveFaceUncoveredFractionsPtr_)
{
calcAddressing();
}
return *slaveFaceUncoveredFractionsPtr_;
}
template<class MasterPatch, class SlavePatch>
bool GGIInterpolation<MasterPatch, SlavePatch>::movePoints
(

49
src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H
View file

@ -219,9 +219,30 @@ class GGIInterpolation
//- List of uncovered master faces
mutable labelList* uncoveredMasterAddrPtr_;
//- List of partially covered master faces. Contains faces with
// weights between masterNonOverlapFaceTol_ and
// uncoveredFaceAreaTol_. Used in bridging
mutable labelList* partiallyCoveredMasterAddrPtr_;
//- List of face fractions for partially covered master faces. Range
// anywhere from 1 to uncoveredFaceAreaTol_ (e.g. 0.2 means
// that the face area is covered by 80% on the other side)
mutable scalarField* masterFaceUncoveredFractionsPtr_;
//- List of uncovered slave faces
mutable labelList* uncoveredSlaveAddrPtr_;
//- List of partially covered slave faces. Contains faces with
// weights between masterNonOverlapFaceTol_ and
// uncoveredFaceAreaTol_. Used in bridging
mutable labelList* partiallyCoveredSlaveAddrPtr_;
//- List of face fractions for partially covered slave faces. Range
// anywhere from 1 to uncoveredFaceAreaTol_ (e.g. 0.2 means
// that the face area is covered by 80% on the other side)
mutable scalarField* slaveFaceUncoveredFractionsPtr_;
// Private static data
//- Facet area error tolerance
@ -233,6 +254,9 @@ class GGIInterpolation
//- Facet bound box extension factor
static const debug::tolerancesSwitch faceBoundBoxExtendSpanFraction_;
//- Tolerance for uncovered face areas
static const debug::tolerancesSwitch uncoveredFaceAreaTol_;
//- The next 3 attributes are parameters controlling the creation
// of an octree search engine for the GGI facets neighbourhood
// determination.
@ -417,6 +441,14 @@ class GGIInterpolation
const scalar& nonOverlapFaceTol
) const;
//- Calculate partially covered faces
void calcPartiallyCoveredFaces
(
const scalarListList& patchWeights,
const scalar& nonOverlapFaceTo,
const bool isMaster
) const;
//- Clear all geometry and addressing
void clearOut();
@ -451,7 +483,9 @@ class GGIInterpolation
(
const Field<Type>& bridgeField,
Field<Type>& ff,
const labelList& addr
const labelList& addr,
const labelList& partiallyCoveredAddr,
const scalarField& uncoveredFractions
);
//- Bridge uncovered faces given addressing
@ -533,9 +567,21 @@ public:
//- Return reference to the master list of non-overlap faces
const labelList& uncoveredMasterFaces() const;
//- Return reference to master list of partially covered faces
const labelList& partiallyCoveredMasterFaces() const;
//- Return uncovered fractions of partially covered master faces
const scalarField& masterFaceUncoveredFractions() const;
//- Return reference to the slave list of non-overlap faces
const labelList& uncoveredSlaveFaces() const;
//- Return reference to slave list of partially covered faces
const labelList& partiallyCoveredSlaveFaces() const;
//- Return uncovered fractions of partially covered slave faces
const scalarField& slaveFaceUncoveredFractions() const;
//- Return reference to point addressing
const List<labelPair>& masterPointAddr() const;
@ -554,6 +600,7 @@ public:
//- Return distance to intersection for patch points
const scalarField& slavePointDistanceToIntersection() const;
// Interpolation functions
//- Interpolate from master to slave

147
src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C
View file

@ -68,6 +68,16 @@ GGIInterpolation<MasterPatch, SlavePatch>::featureCosTol_
);
template<class MasterPatch, class SlavePatch>
const Foam::debug::tolerancesSwitch
GGIInterpolation<MasterPatch, SlavePatch>::uncoveredFaceAreaTol_
(
"GGIUncoveredFaceAreaTol",
0.95,
"Fraction of face area mismatch (sum of weights) to consider a face "
"as uncovered, i.e. not to rescale weights."
);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class MasterPatch, class SlavePatch>
@ -80,7 +90,11 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
|| slaveAddrPtr_
|| slaveWeightsPtr_
|| uncoveredMasterAddrPtr_
|| partiallyCoveredMasterAddrPtr_
|| masterFaceUncoveredFractionsPtr_
|| uncoveredSlaveAddrPtr_
|| partiallyCoveredSlaveAddrPtr_
|| slaveFaceUncoveredFractionsPtr_
)
{
FatalErrorIn
@ -601,6 +615,29 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
findNonOverlappingFaces(saW, slaveNonOverlapFaceTol_)
);
// Calculate master and slave partially covered addressing
// Note: This function allocates:
// 1. partiallyCoveredMasterAddrPtr_
// 2. masterFaceCoveredFractionsPtr_
calcPartiallyCoveredFaces
(
maW,
masterNonOverlapFaceTol_,
true // This is master
);
// Note: this function allocates:
// 1. partiallyCoveredSlaveAddrPtr_
// 2. slaveFaceCoveredFractionsPtr_
calcPartiallyCoveredFaces
(
saW,
slaveNonOverlapFaceTol_,
false // This is not master
);
// Rescaling the weighting factors so they will sum up to 1.0
// See the comment for the method ::rescaleWeightingFactors() for
// more information. By default, we always rescale. But for some
@ -608,7 +645,10 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
// then we need the brute values, so no rescaling in that
// case. Hence the little flag rescaleGGIWeightingFactors_
// Not parallelised. HJ, 27/Apr/2016
// Not parallelised. HJ, 27/Apr/2016. Rescaling is currently switched off
// if the bridge overlap option is used. We should maybe perform rescaling
// for faces that are considered to fully overlap according to
// uncoveredFaceAreaTol_? Reconsider, VV, 12/Oct/2017
if (rescaleGGIWeightingFactors_)
{
rescaleWeightingFactors();
@ -767,6 +807,111 @@ GGIInterpolation<MasterPatch, SlavePatch>::findNonOverlappingFaces
return tpatchFaceNonOverlapAddr;
}
template<class MasterPatch, class SlavePatch>
void GGIInterpolation<MasterPatch, SlavePatch>::calcPartiallyCoveredFaces
(
const scalarListList& patchWeights,
const scalar& nonOverlapFaceTol,
const bool isMaster
) const
{
// Sanity checks first
if
(
isMaster
&& (partiallyCoveredMasterAddrPtr_ || masterFaceUncoveredFractionsPtr_)
)
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"calcPartiallyCoveredFaces() const"
) << "Already calculated master partially covered faces"
<< abort(FatalError);
}
if
(
!isMaster
&& (partiallyCoveredSlaveAddrPtr_ || slaveFaceUncoveredFractionsPtr_)
)
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"calcPartiallyCoveredFaces() const"
) << "Already calculated slave partially covered faces"
<< abort(FatalError);
}
// Temporary storage
DynamicList<label, 64> patchFacePartialOverlap(patchWeights.size());
DynamicList<scalar, 64> uncoveredFaceFractions(patchWeights.size());
// Scan the list of patch weights and collect ones inbetween
// nonOverlapFaceTol and uncoveredFaceAreaTol_
forAll (patchWeights, paWi)
{
const scalar sumWeightsFace = sum(patchWeights[paWi]);
if
(
sumWeightsFace > nonOverlapFaceTol
&& sumWeightsFace <= uncoveredFaceAreaTol_()
)
{
// This is considered partially overlapped face, store the index and
// the non-overlapping area (1 - sum of weights)
patchFacePartialOverlap.append(paWi);
uncoveredFaceFractions.append(1.0 - sumWeightsFace);
}
}
// Transfer the storage
if (isMaster)
{
// Allocate master side
partiallyCoveredMasterAddrPtr_ =
new labelList(patchFacePartialOverlap.xfer());
masterFaceUncoveredFractionsPtr_ =
new scalarField(uncoveredFaceFractions.xfer());
if (debug)
{
InfoIn("GGIInterpolation::calcPartiallyCoveredFaces")
<< " : Found " << partiallyCoveredMasterAddrPtr_->size()
<< " partially overlapping faces for this GGI patch" << nl
<< "Max uncoverage: "
<< max(*masterFaceUncoveredFractionsPtr_)
<< ", min uncoverage: "
<< min(*masterFaceUncoveredFractionsPtr_)
<< endl;
}
}
else
{
// Allocate slave side
partiallyCoveredSlaveAddrPtr_ =
new labelList(patchFacePartialOverlap.xfer());
slaveFaceUncoveredFractionsPtr_ =
new scalarField(uncoveredFaceFractions.xfer());
if (debug)
{
InfoIn("GGIInterpolation::calcPartiallyCoveredFaces")
<< " : Found " << partiallyCoveredSlaveAddrPtr_->size()
<< " partially overlapping faces for this GGI patch" << nl
<< "Max uncoverage: "
<< max(*slaveFaceUncoveredFractionsPtr_)
<< ", min uncoverage: "
<< min(*slaveFaceUncoveredFractionsPtr_)
<< endl;
}
}
}
template<class FromPatch, class ToPatch>
void GGIInterpolation<FromPatch, ToPatch>::
calcMasterPointAddressing() const

6
src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C
View file

@ -227,7 +227,11 @@ void Foam::ggiPolyPatch::calcPatchToPatch() const
// HJ, 30/Jun/2013
SMALL, // Non-overlapping face tolerances
SMALL, // HJ, 24/Oct/2008
true, // Rescale weighting factors. Bug fix, MB.
bridgeOverlap_ ? false : true, // Do not rescale weights if
// bridge overlap is switched on
// as this will mess up
// partially overlapped
// faces. VV, 12/Oct/2017.
reject_ // Quick rejection algorithm, default BB_OCTREE
);