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

Bugfix: partial overlap ggi conservation fix for significantly uncovered partial overlap faces. Vuko Vukcevic

Browse source
This commit is contained in:
Hrvoje Jasak 2017-10-20 10:21:53 +01:00
parent 246a172c9d
commit 255a5c7940
7 changed files with 450 additions and 55 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -163,8 +163,7 @@ tmp<Field<Type> > ggiFvPatchField<Type>::patchNeighbourField() const
// Get shadow face-cells and assemble shadow field
// This is a patchInternalField of neighbour but access is inconvenient.
// Assemble by hand.
// HJ, 27/Sep/2011
// Assemble by hand. HJ, 27/Sep/2011
const unallocLabelList& sfc = ggiPatch_.shadow().faceCells();
Field<Type> sField(sfc.size());
@ -180,13 +179,15 @@ tmp<Field<Type> > ggiFvPatchField<Type>::patchNeighbourField() const
if (ggiPatch_.bridgeOverlap())
{
// Symmetry treatment used for overlap
vectorField nHat = this->patch().nf();
const vectorField nHat = this->patch().nf();
// Use mirrored neighbour field for interpolation
// HJ, 21/Jan/2009
Field<Type> bridgeField =
const Field<Type> bridgeField =
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);
}
@ -211,18 +212,8 @@ void ggiFvPatchField<Type>::initEvaluate
+ (1.0 - this->patch().weights())*this->patchNeighbourField()
);
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);
}
// Note: bridging and correction of partially overlapping faces taken into
// account in patchNeighbourField(). VV, 16/Oct/2017.
Field<Type>::operator=(pf);
}
@ -266,6 +257,21 @@ void ggiFvPatchField<Type>::initInterfaceMatrixUpdate
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
const unallocLabelList& fc = ggiPatch_.faceCells();
@ -327,6 +333,21 @@ void ggiFvPatchField<Type>::initInterfaceMatrixUpdate
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
// of multiplication
multiply(pnf, coeffs, pnf);

44
src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C
View file

@ -62,14 +62,14 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const
// HJ, 2/Aug/2007
if (ggiPolyPatch_.master())
{
vectorField n = nf();
const vectorField n = nf();
// Note: mag in the dot-product.
// For all valid meshes, the non-orthogonality will be less than
// 90 deg and the dot-product will be positive. For invalid
// meshes (d & s <= 0), this will stabilise the calculation
// but the result will be poor. HJ, 24/Aug/2011
scalarField nfc =
const scalarField nfc =
mag(n & (ggiPolyPatch_.reconFaceCellCentres() - Cf()));
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
// for interpolation. HJ, 21/Jan/2009
bridge(scalarField(size(), 0.5), w);
const scalarField bridgedField(size(), 0.5);
bridge(bridgedField, w);
}
}
else
@ -87,7 +89,7 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const
scalarField masterWeights(shadow().size());
shadow().makeWeights(masterWeights);
scalarField oneMinusW = 1 - masterWeights;
const scalarField oneMinusW = 1 - masterWeights;
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
// 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())
{
// 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));
if (bridgeOverlap())
{
scalarField bridgeDeltas = nf() & fvPatch::delta();
bridge(bridgeDeltas, dc);
}
// Note: no need to bridge the overlap since delta already takes it into
// account. VV, 18/Oct/2017.
}
else
{
@ -126,7 +126,11 @@ void Foam::ggiFvPatch::makeDeltaCoeffs(scalarField& dc) const
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);
}
@ -143,8 +147,8 @@ void Foam::ggiFvPatch::makeCorrVecs(vectorField& cv) const
// Calculate correction vectors on coupled patches
const scalarField& patchDeltaCoeffs = deltaCoeffs();
vectorField patchDeltas = delta();
vectorField n = nf();
const vectorField patchDeltas = delta();
const vectorField n = nf();
// If non-orthogonality is over 90 deg, kill correction vector
// HJ, 6/Jan/2011
@ -161,7 +165,10 @@ Foam::tmp<Foam::vectorField> Foam::ggiFvPatch::delta() const
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());
}
@ -177,7 +184,10 @@ Foam::tmp<Foam::vectorField> Foam::ggiFvPatch::delta() const
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());
}

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

@ -97,13 +97,26 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridge
(
const Field<Type>& bridgeField,
Field<Type>& ff,
const labelList& addr
const labelList& addr,
const labelList& partiallyCoveredAddr,
const scalarField& coveredFractions
)
{
// Fully uncovered faces
forAll (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,
Field<Type>& ff,
const labelList& mask,
const labelList& uncoveredFaces
const labelList& uncoveredFaces,
const labelList& partiallyCoveredAddr,
const scalarField& coveredFractions
)
{
// Note: tricky algorithm
@ -131,7 +146,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridge
const label faceI = uncoveredFaces[uncoI];
// Search through the mask
for (; maskAddrI < mask.size(); maskAddrI++)
for (; maskAddrI < mask.size(); ++maskAddrI)
{
if (faceI == mask[maskAddrI])
{
@ -147,7 +162,38 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridge
// Go one back and check for next uncovered face
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;
@ -511,7 +557,8 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridgeMaster
if
(
bridgeField.size() != masterPatch_.size()
|| ff.size() != masterPatch_.size())
|| ff.size() != masterPatch_.size()
)
{
FatalErrorIn
(
@ -527,7 +574,14 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridgeMaster
<< 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"
" const labelList& mask\n"
") const"
) << "given field does not correspond to patch. Patch size: "
) << "given field does not correspond to patch. Patch (mask) size: "
<< masterPatch_.size()
<< " bridge field size: " << bridgeField.size()
<< " field size: " << ff.size()
@ -559,7 +613,15 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeMaster
<< 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);
}
bridge(bridgeField, ff, uncoveredSlaveFaces());
bridge
(
bridgeField,
ff,
uncoveredSlaveFaces(),
partiallyCoveredSlaveFaces(),
slaveFaceCoveredFractions()
);
}
@ -616,7 +685,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeSlave
" Field<Type>& ff\n,"
" const labelList& mask\n"
") const"
) << "given field does not correspond to patch. Patch size: "
) << "given field does not correspond to patch. Patch (mask) size: "
<< slavePatch_.size()
<< " bridge field size: " << bridgeField.size()
<< " field size: " << ff.size()
@ -624,12 +693,18 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeSlave
<< abort(FatalError);
}
maskedBridge(bridgeField, ff, mask, uncoveredSlaveFaces());
maskedBridge
(
bridgeField,
ff,
mask,
uncoveredSlaveFaces(),
partiallyCoveredSlaveFaces(),
slaveFaceCoveredFractions()
);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

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(masterFaceCoveredFractionsPtr_);
deleteDemandDrivenData(uncoveredSlaveAddrPtr_);
deleteDemandDrivenData(partiallyCoveredSlaveAddrPtr_);
deleteDemandDrivenData(slaveFaceCoveredFractionsPtr_);
deleteDemandDrivenData(masterPointAddressingPtr_);
deleteDemandDrivenData(masterPointWeightsPtr_);
@ -155,7 +159,11 @@ GGIInterpolation<MasterPatch, SlavePatch>::GGIInterpolation
slavePointWeightsPtr_(NULL),
slavePointDistancePtr_(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
// 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>
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>::slaveFaceCoveredFractions() const
{
if (!slaveFaceCoveredFractionsPtr_)
{
calcAddressing();
}
return *slaveFaceCoveredFractionsPtr_;
}
template<class MasterPatch, class SlavePatch>
bool GGIInterpolation<MasterPatch, SlavePatch>::movePoints
(

53
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 20% on the other side)
mutable scalarField* masterFaceCoveredFractionsPtr_;
//- 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 20% on the other side)
mutable scalarField* slaveFaceCoveredFractionsPtr_;
// 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& coveredFractions
);
//- Bridge uncovered faces given addressing
@ -462,7 +496,9 @@ class GGIInterpolation
const Field<Type>& bridgeField,
Field<Type>& ff,
const labelList& mask,
const labelList& uncoveredFaces
const labelList& uncoveredFaces,
const labelList& partiallyCoveredAddr,
const scalarField& coveredFractions
);
//- Is a transform required?
@ -533,9 +569,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 covered fractions of partially covered master faces
const scalarField& masterFaceCoveredFractions() 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 covered fractions of partially covered slave faces
const scalarField& slaveFaceCoveredFractions() const;
//- Return reference to point addressing
const List<labelPair>& masterPointAddr() const;
@ -554,6 +602,7 @@ public:
//- Return distance to intersection for patch points
const scalarField& slavePointDistanceToIntersection() const;
// Interpolation functions
//- Interpolate from master to slave

190
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.999,
"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_
|| masterFaceCoveredFractionsPtr_
|| uncoveredSlaveAddrPtr_
|| partiallyCoveredSlaveAddrPtr_
|| slaveFaceCoveredFractionsPtr_
)
{
FatalErrorIn
@ -431,7 +445,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
);
// Fix: previously checked for VSMALL.
// HJ, 19/Sep2/106
// HJ, 19/Sep/2016
if (intersectionArea > intersectionTestArea)
{
// We compute the GGI weights based on this
@ -555,7 +569,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
const labelList& curMa = ma[mfI];
const scalarList& cursmaW = smaW[mfI];
// Current master face indes
// Current master face index
const label faceMaster = mfI;
forAll (curMa, mAI)
@ -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,8 @@ 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 not performed for
// partially overlapping faces for their correct treatment. VV, 16/Oct/2017.
if (rescaleGGIWeightingFactors_)
{
rescaleWeightingFactors();
@ -656,6 +694,35 @@ void GGIInterpolation<MasterPatch, SlavePatch>::rescaleWeightingFactors() const
scalar sumMWC = 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
if (debug)
{
@ -679,7 +746,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::rescaleWeightingFactors() const
{
scalar slaveWeightSum = Foam::sum(saW[saWi]);
if (saW[saWi].size() > 0)
if (saW[saWi].size() > 0 && !slavePCMask[saWi])
{
saW[saWi] = saW[saWi]/slaveWeightSum;
@ -696,7 +763,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::rescaleWeightingFactors() const
{
scalar masterWeightSum = Foam::sum(maW[maWi]);
if (maW[maWi].size() > 0)
if (maW[maWi].size() > 0 && !masterPCMask[maWi])
{
maW[maWi] = maW[maWi]/masterWeightSum;
@ -767,6 +834,119 @@ 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_ || 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>
void GGIInterpolation<FromPatch, ToPatch>::
calcMasterPointAddressing() const

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

@ -324,7 +324,7 @@ void Foam::ggiPolyPatch::calcLocalParallel() const
{
FatalErrorIn("void ggiPolyPatch::calcLocalParallel() const")
<< "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 "
<< "global in parallel runs"
<< abort(FatalError);