Merge branch 'bugfix/partialOveralpGGI'

This commit is contained in:
Hrvoje Jasak 2017-10-20 10:22:12 +01:00
commit f500917045
7 changed files with 450 additions and 55 deletions

View file

@ -163,8 +163,7 @@ tmp<Field<Type> > ggiFvPatchField<Type>::patchNeighbourField() const
// Get shadow face-cells and assemble shadow field // Get shadow face-cells and assemble shadow field
// This is a patchInternalField of neighbour but access is inconvenient. // This is a patchInternalField of neighbour but access is inconvenient.
// Assemble by hand. // Assemble by hand. HJ, 27/Sep/2011
// HJ, 27/Sep/2011
const unallocLabelList& sfc = ggiPatch_.shadow().faceCells(); const unallocLabelList& sfc = ggiPatch_.shadow().faceCells();
Field<Type> sField(sfc.size()); Field<Type> sField(sfc.size());
@ -180,13 +179,15 @@ tmp<Field<Type> > ggiFvPatchField<Type>::patchNeighbourField() const
if (ggiPatch_.bridgeOverlap()) if (ggiPatch_.bridgeOverlap())
{ {
// Symmetry treatment used for overlap // Symmetry treatment used for overlap
vectorField nHat = this->patch().nf(); const vectorField nHat = this->patch().nf();
// Use mirrored neighbour field for interpolation // Use mirrored neighbour field for interpolation
// HJ, 21/Jan/2009 // HJ, 21/Jan/2009
Field<Type> bridgeField = const Field<Type> bridgeField =
transform(I - 2.0*sqr(nHat), this->patchInternalField()); transform(I - 2.0*sqr(nHat), this->patchInternalField());
// Note: bridging now takes into account fully uncovered and partially
// covered faces. VV, 18/Oct/2017.
ggiPatch_.bridge(bridgeField, pnf); ggiPatch_.bridge(bridgeField, pnf);
} }
@ -211,18 +212,8 @@ void ggiFvPatchField<Type>::initEvaluate
+ (1.0 - this->patch().weights())*this->patchNeighbourField() + (1.0 - this->patch().weights())*this->patchNeighbourField()
); );
if (ggiPatch_.bridgeOverlap()) // Note: bridging and correction of partially overlapping faces taken into
{ // account in patchNeighbourField(). VV, 16/Oct/2017.
// 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); Field<Type>::operator=(pf);
} }
@ -266,6 +257,21 @@ void ggiFvPatchField<Type>::initInterfaceMatrixUpdate
scalarField pnf = ggiPatch_.interpolate(sField); scalarField pnf = ggiPatch_.interpolate(sField);
if (ggiPatch_.bridgeOverlap())
{
// Note: will not work properly for types with rank > 0 (everything
// above scalar) if the symmetry plane is not aligned with one of the
// coordinate axes. VV, 18/Oct/2017.
const scalarField bridgeField =
transform
(
I - 2.0*sqr(this->patch().nf()),
ggiPatch_.patchInternalField(psiInternal)
);
ggiPatch_.bridge(bridgeField, pnf);
}
// Multiply the field by coefficients and add into the result // Multiply the field by coefficients and add into the result
const unallocLabelList& fc = ggiPatch_.faceCells(); const unallocLabelList& fc = ggiPatch_.faceCells();
@ -327,6 +333,21 @@ void ggiFvPatchField<Type>::initInterfaceMatrixUpdate
Field<Type> pnf = ggiPatch_.interpolate(sField); Field<Type> pnf = ggiPatch_.interpolate(sField);
if (ggiPatch_.bridgeOverlap())
{
// Note: will not work properly for types with rank > 0 (everything
// above scalar) if the symmetry plane is not aligned with one of the
// coordinate axes. VV, 18/Oct/2017.
const Field<Type> bridgeField =
transform
(
I - 2.0*sqr(this->patch().nf()),
ggiPatch_.patchInternalField(psiInternal)
);
ggiPatch_.bridge(bridgeField, pnf);
}
// Multiply neighbour field with coeffs and re-use pnf for result // Multiply neighbour field with coeffs and re-use pnf for result
// of multiplication // of multiplication
multiply(pnf, coeffs, pnf); multiply(pnf, coeffs, pnf);

View file

@ -62,14 +62,14 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const
// HJ, 2/Aug/2007 // HJ, 2/Aug/2007
if (ggiPolyPatch_.master()) if (ggiPolyPatch_.master())
{ {
vectorField n = nf(); const vectorField n = nf();
// Note: mag in the dot-product. // Note: mag in the dot-product.
// For all valid meshes, the non-orthogonality will be less than // For all valid meshes, the non-orthogonality will be less than
// 90 deg and the dot-product will be positive. For invalid // 90 deg and the dot-product will be positive. For invalid
// meshes (d & s <= 0), this will stabilise the calculation // meshes (d & s <= 0), this will stabilise the calculation
// but the result will be poor. HJ, 24/Aug/2011 // but the result will be poor. HJ, 24/Aug/2011
scalarField nfc = const scalarField nfc =
mag(n & (ggiPolyPatch_.reconFaceCellCentres() - Cf())); mag(n & (ggiPolyPatch_.reconFaceCellCentres() - Cf()));
w = nfc/(mag(n & (Cf() - Cn())) + nfc); w = nfc/(mag(n & (Cf() - Cn())) + nfc);
@ -78,7 +78,9 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const
{ {
// Set overlap weights to 0.5 and use mirrored neighbour field // Set overlap weights to 0.5 and use mirrored neighbour field
// for interpolation. HJ, 21/Jan/2009 // for interpolation. HJ, 21/Jan/2009
bridge(scalarField(size(), 0.5), w); const scalarField bridgedField(size(), 0.5);
bridge(bridgedField, w);
} }
} }
else else
@ -87,7 +89,7 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const
scalarField masterWeights(shadow().size()); scalarField masterWeights(shadow().size());
shadow().makeWeights(masterWeights); shadow().makeWeights(masterWeights);
scalarField oneMinusW = 1 - masterWeights; const scalarField oneMinusW = 1 - masterWeights;
w = interpolate(oneMinusW); w = interpolate(oneMinusW);
@ -95,7 +97,9 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const
{ {
// Set overlap weights to 0.5 and use mirrored neighbour field // Set overlap weights to 0.5 and use mirrored neighbour field
// for interpolation. HJ, 21/Jan/2009 // for interpolation. HJ, 21/Jan/2009
bridge(scalarField(size(), 0.5), w); const scalarField bridgedField(size(), 0.5);
bridge(bridgedField, w);
} }
} }
} }
@ -107,16 +111,12 @@ void Foam::ggiFvPatch::makeDeltaCoeffs(scalarField& dc) const
if (ggiPolyPatch_.master()) if (ggiPolyPatch_.master())
{ {
// Stabilised form for bad meshes. HJ, 24/Aug/2011 // Stabilised form for bad meshes. HJ, 24/Aug/2011
vectorField d = delta(); const vectorField d = delta();
dc = 1.0/max(nf() & d, 0.05*mag(d)); dc = 1.0/max(nf() & d, 0.05*mag(d));
if (bridgeOverlap()) // Note: no need to bridge the overlap since delta already takes it into
{ // account. VV, 18/Oct/2017.
scalarField bridgeDeltas = nf() & fvPatch::delta();
bridge(bridgeDeltas, dc);
}
} }
else else
{ {
@ -126,7 +126,11 @@ void Foam::ggiFvPatch::makeDeltaCoeffs(scalarField& dc) const
if (bridgeOverlap()) if (bridgeOverlap())
{ {
scalarField bridgeDeltas = nf() & fvPatch::delta(); // Note: double the deltaCoeffs because this is symmetry treatment
// and fvPatch::deltaCoeffs() is cell to face inverse distance,
// while we need cell to "symmetry neighbour cell" distance.
// VV, 18/Oct/2017.
const scalarField bridgeDeltas = 2.0*fvPatch::deltaCoeffs();
bridge(bridgeDeltas, dc); bridge(bridgeDeltas, dc);
} }
@ -143,8 +147,8 @@ void Foam::ggiFvPatch::makeCorrVecs(vectorField& cv) const
// Calculate correction vectors on coupled patches // Calculate correction vectors on coupled patches
const scalarField& patchDeltaCoeffs = deltaCoeffs(); const scalarField& patchDeltaCoeffs = deltaCoeffs();
vectorField patchDeltas = delta(); const vectorField patchDeltas = delta();
vectorField n = nf(); const vectorField n = nf();
// If non-orthogonality is over 90 deg, kill correction vector // If non-orthogonality is over 90 deg, kill correction vector
// HJ, 6/Jan/2011 // HJ, 6/Jan/2011
@ -161,7 +165,10 @@ Foam::tmp<Foam::vectorField> Foam::ggiFvPatch::delta() const
if (bridgeOverlap()) if (bridgeOverlap())
{ {
vectorField bridgeDeltas = Cf() - Cn(); // Note: double the deltas because this is symmetry treatment and
// fvPatch::delta() is cell to face distance, while we need cell to
// "symmetry neighbour cell" distance. VV, 18/Oct/2017.
const vectorField bridgeDeltas = 2.0*fvPatch::delta();
bridge(bridgeDeltas, tDelta()); bridge(bridgeDeltas, tDelta());
} }
@ -177,7 +184,10 @@ Foam::tmp<Foam::vectorField> Foam::ggiFvPatch::delta() const
if (bridgeOverlap()) if (bridgeOverlap())
{ {
vectorField bridgeDeltas = Cf() - Cn(); // Note: double the deltas because this is symmetry treatment and
// fvPatch::delta() is cell to face distance, while we need cell to
// "symmetry neighbour cell" distance. VV, 18/Oct/2017.
const vectorField bridgeDeltas = 2.0*fvPatch::delta();
bridge(bridgeDeltas, tDelta()); bridge(bridgeDeltas, tDelta());
} }

View file

@ -97,13 +97,26 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridge
( (
const Field<Type>& bridgeField, const Field<Type>& bridgeField,
Field<Type>& ff, Field<Type>& ff,
const labelList& addr const labelList& addr,
const labelList& partiallyCoveredAddr,
const scalarField& coveredFractions
) )
{ {
// Fully uncovered faces
forAll (addr, faceI) forAll (addr, faceI)
{ {
ff[addr[faceI]] = bridgeField[addr[faceI]]; ff[addr[faceI]] = bridgeField[addr[faceI]];
} }
// Loop through partially covered faces and correct them. Note the
// operator+= since we assume that the interpolation part is carried out
// before bridging (see e.g. ggiFvPatchField::patchNeighbourField()) using
// weights that do not sum up to 1
forAll (partiallyCoveredAddr, pcfI)
{
ff[partiallyCoveredAddr[pcfI]] +=
coveredFractions[pcfI]*bridgeField[partiallyCoveredAddr[pcfI]];
}
} }
@ -114,7 +127,9 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridge
const Field<Type>& bridgeField, const Field<Type>& bridgeField,
Field<Type>& ff, Field<Type>& ff,
const labelList& mask, const labelList& mask,
const labelList& uncoveredFaces const labelList& uncoveredFaces,
const labelList& partiallyCoveredAddr,
const scalarField& coveredFractions
) )
{ {
// Note: tricky algorithm // Note: tricky algorithm
@ -131,7 +146,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridge
const label faceI = uncoveredFaces[uncoI]; const label faceI = uncoveredFaces[uncoI];
// Search through the mask // Search through the mask
for (; maskAddrI < mask.size(); maskAddrI++) for (; maskAddrI < mask.size(); ++maskAddrI)
{ {
if (faceI == mask[maskAddrI]) if (faceI == mask[maskAddrI])
{ {
@ -147,7 +162,38 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridge
// Go one back and check for next uncovered face // Go one back and check for next uncovered face
if (maskAddrI > 0) if (maskAddrI > 0)
{ {
maskAddrI--; --maskAddrI;
}
break;
}
}
}
// Reset maskAddrI
maskAddrI = 0;
forAll (partiallyCoveredAddr, pcfI)
{
// Pick partially covered face
const label faceI = partiallyCoveredAddr[pcfI];
for (; maskAddrI < mask.size(); ++maskAddrI)
{
if (faceI == mask[maskAddrI])
{
// Found masked partially covered face
ff[maskAddrI] += coveredFractions[pcfI]*bridgeField[maskAddrI];
break;
}
else if (mask[maskAddrI] > faceI)
{
// Gone beyond my index: my face is not present in the mask
// Go one back and check for next uncovered face
if (maskAddrI > 0)
{
--maskAddrI;
} }
break; break;
@ -511,7 +557,8 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridgeMaster
if if
( (
bridgeField.size() != masterPatch_.size() bridgeField.size() != masterPatch_.size()
|| ff.size() != masterPatch_.size()) || ff.size() != masterPatch_.size()
)
{ {
FatalErrorIn FatalErrorIn
( (
@ -527,7 +574,14 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridgeMaster
<< abort(FatalError); << abort(FatalError);
} }
bridge(bridgeField, ff, uncoveredMasterFaces()); bridge
(
bridgeField,
ff,
uncoveredMasterFaces(),
partiallyCoveredMasterFaces(),
masterFaceCoveredFractions()
);
} }
@ -551,7 +605,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeMaster
" Field<Type>& ff,\n" " Field<Type>& ff,\n"
" const labelList& mask\n" " const labelList& mask\n"
") const" ") const"
) << "given field does not correspond to patch. Patch size: " ) << "given field does not correspond to patch. Patch (mask) size: "
<< masterPatch_.size() << masterPatch_.size()
<< " bridge field size: " << bridgeField.size() << " bridge field size: " << bridgeField.size()
<< " field size: " << ff.size() << " field size: " << ff.size()
@ -559,7 +613,15 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeMaster
<< abort(FatalError); << abort(FatalError);
} }
maskedBridge(bridgeField, ff, mask, uncoveredMasterFaces()); maskedBridge
(
bridgeField,
ff,
mask,
uncoveredMasterFaces(),
partiallyCoveredMasterFaces(),
masterFaceCoveredFractions()
);
} }
@ -592,7 +654,14 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridgeSlave
<< abort(FatalError); << abort(FatalError);
} }
bridge(bridgeField, ff, uncoveredSlaveFaces()); bridge
(
bridgeField,
ff,
uncoveredSlaveFaces(),
partiallyCoveredSlaveFaces(),
slaveFaceCoveredFractions()
);
} }
@ -616,7 +685,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeSlave
" Field<Type>& ff\n," " Field<Type>& ff\n,"
" const labelList& mask\n" " const labelList& mask\n"
") const" ") const"
) << "given field does not correspond to patch. Patch size: " ) << "given field does not correspond to patch. Patch (mask) size: "
<< slavePatch_.size() << slavePatch_.size()
<< " bridge field size: " << bridgeField.size() << " bridge field size: " << bridgeField.size()
<< " field size: " << ff.size() << " field size: " << ff.size()
@ -624,12 +693,18 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeSlave
<< abort(FatalError); << abort(FatalError);
} }
maskedBridge(bridgeField, ff, mask, uncoveredSlaveFaces()); maskedBridge
(
bridgeField,
ff,
mask,
uncoveredSlaveFaces(),
partiallyCoveredSlaveFaces(),
slaveFaceCoveredFractions()
);
} }
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam } // End namespace Foam

View file

@ -104,7 +104,11 @@ void GGIInterpolation<MasterPatch, SlavePatch>::clearOut()
deleteDemandDrivenData(slaveWeightsPtr_); deleteDemandDrivenData(slaveWeightsPtr_);
deleteDemandDrivenData(uncoveredMasterAddrPtr_); deleteDemandDrivenData(uncoveredMasterAddrPtr_);
deleteDemandDrivenData(partiallyCoveredMasterAddrPtr_);
deleteDemandDrivenData(masterFaceCoveredFractionsPtr_);
deleteDemandDrivenData(uncoveredSlaveAddrPtr_); deleteDemandDrivenData(uncoveredSlaveAddrPtr_);
deleteDemandDrivenData(partiallyCoveredSlaveAddrPtr_);
deleteDemandDrivenData(slaveFaceCoveredFractionsPtr_);
deleteDemandDrivenData(masterPointAddressingPtr_); deleteDemandDrivenData(masterPointAddressingPtr_);
deleteDemandDrivenData(masterPointWeightsPtr_); deleteDemandDrivenData(masterPointWeightsPtr_);
@ -155,7 +159,11 @@ GGIInterpolation<MasterPatch, SlavePatch>::GGIInterpolation
slavePointWeightsPtr_(NULL), slavePointWeightsPtr_(NULL),
slavePointDistancePtr_(NULL), slavePointDistancePtr_(NULL),
uncoveredMasterAddrPtr_(NULL), uncoveredMasterAddrPtr_(NULL),
uncoveredSlaveAddrPtr_(NULL) partiallyCoveredMasterAddrPtr_(NULL),
masterFaceCoveredFractionsPtr_(NULL),
uncoveredSlaveAddrPtr_(NULL),
partiallyCoveredSlaveAddrPtr_(NULL),
slaveFaceCoveredFractionsPtr_(NULL)
{ {
// Check size of transform. They should be equal to slave patch size // Check size of transform. They should be equal to slave patch size
// if the transform is not constant // 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>::masterFaceCoveredFractions() const
{
if (!masterFaceCoveredFractionsPtr_)
{
calcAddressing();
}
return *masterFaceCoveredFractionsPtr_;
}
template<class MasterPatch, class SlavePatch> template<class MasterPatch, class SlavePatch>
const labelList& const labelList&
GGIInterpolation<MasterPatch, SlavePatch>::uncoveredSlaveFaces() const 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>::slaveFaceCoveredFractions() const
{
if (!slaveFaceCoveredFractionsPtr_)
{
calcAddressing();
}
return *slaveFaceCoveredFractionsPtr_;
}
template<class MasterPatch, class SlavePatch> template<class MasterPatch, class SlavePatch>
bool GGIInterpolation<MasterPatch, SlavePatch>::movePoints bool GGIInterpolation<MasterPatch, SlavePatch>::movePoints
( (

View file

@ -219,9 +219,30 @@ class GGIInterpolation
//- List of uncovered master faces //- List of uncovered master faces
mutable labelList* uncoveredMasterAddrPtr_; 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 20% on the other side)
mutable scalarField* masterFaceCoveredFractionsPtr_;
//- List of uncovered slave faces //- List of uncovered slave faces
mutable labelList* uncoveredSlaveAddrPtr_; 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 20% on the other side)
mutable scalarField* slaveFaceCoveredFractionsPtr_;
// Private static data // Private static data
//- Facet area error tolerance //- Facet area error tolerance
@ -233,6 +254,9 @@ class GGIInterpolation
//- Facet bound box extension factor //- Facet bound box extension factor
static const debug::tolerancesSwitch faceBoundBoxExtendSpanFraction_; static const debug::tolerancesSwitch faceBoundBoxExtendSpanFraction_;
//- Tolerance for uncovered face areas
static const debug::tolerancesSwitch uncoveredFaceAreaTol_;
//- The next 3 attributes are parameters controlling the creation //- The next 3 attributes are parameters controlling the creation
// of an octree search engine for the GGI facets neighbourhood // of an octree search engine for the GGI facets neighbourhood
// determination. // determination.
@ -417,6 +441,14 @@ class GGIInterpolation
const scalar& nonOverlapFaceTol const scalar& nonOverlapFaceTol
) const; ) const;
//- Calculate partially covered faces
void calcPartiallyCoveredFaces
(
const scalarListList& patchWeights,
const scalar& nonOverlapFaceTo,
const bool isMaster
) const;
//- Clear all geometry and addressing //- Clear all geometry and addressing
void clearOut(); void clearOut();
@ -451,7 +483,9 @@ class GGIInterpolation
( (
const Field<Type>& bridgeField, const Field<Type>& bridgeField,
Field<Type>& ff, Field<Type>& ff,
const labelList& addr const labelList& addr,
const labelList& partiallyCoveredAddr,
const scalarField& coveredFractions
); );
//- Bridge uncovered faces given addressing //- Bridge uncovered faces given addressing
@ -462,7 +496,9 @@ class GGIInterpolation
const Field<Type>& bridgeField, const Field<Type>& bridgeField,
Field<Type>& ff, Field<Type>& ff,
const labelList& mask, const labelList& mask,
const labelList& uncoveredFaces const labelList& uncoveredFaces,
const labelList& partiallyCoveredAddr,
const scalarField& coveredFractions
); );
//- Is a transform required? //- Is a transform required?
@ -533,9 +569,21 @@ public:
//- Return reference to the master list of non-overlap faces //- Return reference to the master list of non-overlap faces
const labelList& uncoveredMasterFaces() const; const labelList& uncoveredMasterFaces() const;
//- Return reference to master list of partially covered faces
const labelList& partiallyCoveredMasterFaces() const;
//- Return covered fractions of partially covered master faces
const scalarField& masterFaceCoveredFractions() const;
//- Return reference to the slave list of non-overlap faces //- Return reference to the slave list of non-overlap faces
const labelList& uncoveredSlaveFaces() const; const labelList& uncoveredSlaveFaces() const;
//- Return reference to slave list of partially covered faces
const labelList& partiallyCoveredSlaveFaces() const;
//- Return covered fractions of partially covered slave faces
const scalarField& slaveFaceCoveredFractions() const;
//- Return reference to point addressing //- Return reference to point addressing
const List<labelPair>& masterPointAddr() const; const List<labelPair>& masterPointAddr() const;
@ -554,6 +602,7 @@ public:
//- Return distance to intersection for patch points //- Return distance to intersection for patch points
const scalarField& slavePointDistanceToIntersection() const; const scalarField& slavePointDistanceToIntersection() const;
// Interpolation functions // Interpolation functions
//- Interpolate from master to slave //- Interpolate from master to slave

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.999,
"Fraction of face area mismatch (sum of weights) to consider a face "
"as uncovered, i.e. not to rescale weights."
);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class MasterPatch, class SlavePatch> template<class MasterPatch, class SlavePatch>
@ -80,7 +90,11 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
|| slaveAddrPtr_ || slaveAddrPtr_
|| slaveWeightsPtr_ || slaveWeightsPtr_
|| uncoveredMasterAddrPtr_ || uncoveredMasterAddrPtr_
|| partiallyCoveredMasterAddrPtr_
|| masterFaceCoveredFractionsPtr_
|| uncoveredSlaveAddrPtr_ || uncoveredSlaveAddrPtr_
|| partiallyCoveredSlaveAddrPtr_
|| slaveFaceCoveredFractionsPtr_
) )
{ {
FatalErrorIn FatalErrorIn
@ -431,7 +445,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
); );
// Fix: previously checked for VSMALL. // Fix: previously checked for VSMALL.
// HJ, 19/Sep2/106 // HJ, 19/Sep/2016
if (intersectionArea > intersectionTestArea) if (intersectionArea > intersectionTestArea)
{ {
// We compute the GGI weights based on this // We compute the GGI weights based on this
@ -555,7 +569,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
const labelList& curMa = ma[mfI]; const labelList& curMa = ma[mfI];
const scalarList& cursmaW = smaW[mfI]; const scalarList& cursmaW = smaW[mfI];
// Current master face indes // Current master face index
const label faceMaster = mfI; const label faceMaster = mfI;
forAll (curMa, mAI) forAll (curMa, mAI)
@ -601,6 +615,29 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
findNonOverlappingFaces(saW, slaveNonOverlapFaceTol_) 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 // Rescaling the weighting factors so they will sum up to 1.0
// See the comment for the method ::rescaleWeightingFactors() for // See the comment for the method ::rescaleWeightingFactors() for
// more information. By default, we always rescale. But for some // more information. By default, we always rescale. But for some
@ -608,7 +645,8 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
// then we need the brute values, so no rescaling in that // then we need the brute values, so no rescaling in that
// case. Hence the little flag rescaleGGIWeightingFactors_ // case. Hence the little flag rescaleGGIWeightingFactors_
// Not parallelised. HJ, 27/Apr/2016 // Not parallelised. HJ, 27/Apr/2016. Rescaling is not performed for
// partially overlapping faces for their correct treatment. VV, 16/Oct/2017.
if (rescaleGGIWeightingFactors_) if (rescaleGGIWeightingFactors_)
{ {
rescaleWeightingFactors(); rescaleWeightingFactors();
@ -656,6 +694,35 @@ void GGIInterpolation<MasterPatch, SlavePatch>::rescaleWeightingFactors() const
scalar sumMWC = 0; scalar sumMWC = 0;
scalar curMWC = 0; scalar curMWC = 0;
// Note: do not rescale weighting factors for partially covered faces
if (!partiallyCoveredMasterAddrPtr_ || !partiallyCoveredSlaveAddrPtr_)
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"rescaleWeightingFactors() const"
) << "Master or slave partially covered faces are not calculated."
<< abort(FatalError);
}
const labelList& partiallyCoveredMasterFaces =
*partiallyCoveredMasterAddrPtr_;
const labelList& partiallyCoveredSlaveFaces =
*partiallyCoveredSlaveAddrPtr_;
// Create a mask for partially covered master/slave faces
boolList masterPCMask(maW.size(), false);
boolList slavePCMask(saW.size(), false);
forAll (partiallyCoveredMasterFaces, pfmI)
{
masterPCMask[partiallyCoveredMasterFaces[pfmI]] = true;
}
forAll (partiallyCoveredSlaveFaces, pfsI)
{
slavePCMask[partiallyCoveredSlaveFaces[pfsI]] = true;
}
// Rescaling the slave weights // Rescaling the slave weights
if (debug) if (debug)
{ {
@ -679,7 +746,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::rescaleWeightingFactors() const
{ {
scalar slaveWeightSum = Foam::sum(saW[saWi]); scalar slaveWeightSum = Foam::sum(saW[saWi]);
if (saW[saWi].size() > 0) if (saW[saWi].size() > 0 && !slavePCMask[saWi])
{ {
saW[saWi] = saW[saWi]/slaveWeightSum; saW[saWi] = saW[saWi]/slaveWeightSum;
@ -696,7 +763,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::rescaleWeightingFactors() const
{ {
scalar masterWeightSum = Foam::sum(maW[maWi]); scalar masterWeightSum = Foam::sum(maW[maWi]);
if (maW[maWi].size() > 0) if (maW[maWi].size() > 0 && !masterPCMask[maWi])
{ {
maW[maWi] = maW[maWi]/masterWeightSum; maW[maWi] = maW[maWi]/masterWeightSum;
@ -767,6 +834,119 @@ GGIInterpolation<MasterPatch, SlavePatch>::findNonOverlappingFaces
return tpatchFaceNonOverlapAddr; 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_ || masterFaceCoveredFractionsPtr_)
)
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"calcPartiallyCoveredFaces() const"
) << "Already calculated master partially covered faces"
<< abort(FatalError);
}
if
(
!isMaster
&& (partiallyCoveredSlaveAddrPtr_ || slaveFaceCoveredFractionsPtr_)
)
{
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());
masterFaceCoveredFractionsPtr_ =
new scalarField(uncoveredFaceFractions.xfer());
if (debug)
{
InfoIn("GGIInterpolation::calcPartiallyCoveredFaces")
<< " : Found " << partiallyCoveredMasterAddrPtr_->size()
<< " partially overlapping faces for master GGI patch" << endl;
if (partiallyCoveredMasterAddrPtr_->size())
{
Info<< "Max coverage: "
<< max(*masterFaceCoveredFractionsPtr_)
<< ", min coverage: "
<< min(*masterFaceCoveredFractionsPtr_)
<< endl;
}
}
}
else
{
// Allocate slave side
partiallyCoveredSlaveAddrPtr_ =
new labelList(patchFacePartialOverlap.xfer());
slaveFaceCoveredFractionsPtr_ =
new scalarField(uncoveredFaceFractions.xfer());
if (debug)
{
InfoIn("GGIInterpolation::calcPartiallyCoveredFaces")
<< " : Found " << partiallyCoveredSlaveAddrPtr_->size()
<< " partially overlapping faces for slave GGI patch" << endl;
if (partiallyCoveredSlaveAddrPtr_->size())
{
Info<< "Max coverage: "
<< max(*slaveFaceCoveredFractionsPtr_)
<< ", min coverage: "
<< min(*slaveFaceCoveredFractionsPtr_)
<< endl;
}
}
}
}
template<class FromPatch, class ToPatch> template<class FromPatch, class ToPatch>
void GGIInterpolation<FromPatch, ToPatch>:: void GGIInterpolation<FromPatch, ToPatch>::
calcMasterPointAddressing() const calcMasterPointAddressing() const

View file

@ -324,7 +324,7 @@ void Foam::ggiPolyPatch::calcLocalParallel() const
{ {
FatalErrorIn("void ggiPolyPatch::calcLocalParallel() const") FatalErrorIn("void ggiPolyPatch::calcLocalParallel() const")
<< "Patch size is greater than zone size for GGI patch " << "Patch size is greater than zone size for GGI patch "
<< name() << ". This is not allowerd: " << name() << ". This is not allowed: "
<< "the face zone must contain all patch faces and be " << "the face zone must contain all patch faces and be "
<< "global in parallel runs" << "global in parallel runs"
<< abort(FatalError); << abort(FatalError);