From ba798d60aca70b1508097348bc464e80da506538 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Fri, 13 Oct 2017 08:35:13 +0200 Subject: [PATCH] Updates to GGIInterpolation Correct handling of partially overlapped faces. Initial commit, not tested. --- .../constraint/ggi/ggiFvPatchField.C | 24 +-- .../GGIInterpolation/GGIInterpolate.C | 31 +++- .../GGIInterpolation/GGIInterpolation.C | 62 +++++++- .../GGIInterpolationTemplate.H | 49 +++++- .../GGIInterpolationWeights.C | 147 +++++++++++++++++- .../polyPatches/constraint/ggi/ggiPolyPatch.C | 6 +- 6 files changed, 301 insertions(+), 18 deletions(-) diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C index fe4628f75..d3e7a11da 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C @@ -211,18 +211,20 @@ void ggiFvPatchField::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 pif = this->patchInternalField(); - - Field 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 pif = this->patchInternalField(); +// +// Field bridgeField = +// 0.5*(pif + transform(I - 2.0*sqr(nHat), pif)); +// +// ggiPatch_.bridge(bridgeField, pf); +// } Field::operator=(pf); } diff --git a/src/foam/interpolations/GGIInterpolation/GGIInterpolate.C b/src/foam/interpolations/GGIInterpolation/GGIInterpolate.C index b5fa94572..8c969b849 100644 --- a/src/foam/interpolations/GGIInterpolation/GGIInterpolate.C +++ b/src/foam/interpolations/GGIInterpolation/GGIInterpolate.C @@ -97,13 +97,24 @@ void GGIInterpolation::bridge ( const Field& bridgeField, Field& 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::bridgeMaster << abort(FatalError); } - bridge(bridgeField, ff, uncoveredMasterFaces()); + bridge + ( + bridgeField, + ff, + uncoveredMasterFaces(), + partiallyCoveredMasterFaces(), + masterFaceUncoveredFractions() + ); } @@ -592,7 +610,14 @@ void GGIInterpolation::bridgeSlave << abort(FatalError); } - bridge(bridgeField, ff, uncoveredSlaveFaces()); + bridge + ( + bridgeField, + ff, + uncoveredSlaveFaces(), + partiallyCoveredSlaveFaces(), + slaveFaceUncoveredFractions() + ); } diff --git a/src/foam/interpolations/GGIInterpolation/GGIInterpolation.C b/src/foam/interpolations/GGIInterpolation/GGIInterpolation.C index 3bf5d55bd..a26b15370 100644 --- a/src/foam/interpolations/GGIInterpolation/GGIInterpolation.C +++ b/src/foam/interpolations/GGIInterpolation/GGIInterpolation.C @@ -104,7 +104,11 @@ void GGIInterpolation::clearOut() deleteDemandDrivenData(slaveWeightsPtr_); deleteDemandDrivenData(uncoveredMasterAddrPtr_); + deleteDemandDrivenData(partiallyCoveredMasterAddrPtr_); + deleteDemandDrivenData(masterFaceUncoveredFractionsPtr_); deleteDemandDrivenData(uncoveredSlaveAddrPtr_); + deleteDemandDrivenData(partiallyCoveredSlaveAddrPtr_); + deleteDemandDrivenData(slaveFaceUncoveredFractionsPtr_); deleteDemandDrivenData(masterPointAddressingPtr_); deleteDemandDrivenData(masterPointWeightsPtr_); @@ -155,7 +159,11 @@ GGIInterpolation::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::uncoveredMasterFaces() const } +template +const labelList& +GGIInterpolation::partiallyCoveredMasterFaces() const +{ + if (!partiallyCoveredMasterAddrPtr_) + { + calcAddressing(); + } + + return *partiallyCoveredMasterAddrPtr_; +} + + +template +const scalarField& +GGIInterpolation::masterFaceUncoveredFractions() const +{ + if (!masterFaceUncoveredFractionsPtr_) + { + calcAddressing(); + } + + return *masterFaceUncoveredFractionsPtr_; +} + + template const labelList& GGIInterpolation::uncoveredSlaveFaces() const @@ -269,6 +303,32 @@ GGIInterpolation::uncoveredSlaveFaces() const } +template +const labelList& +GGIInterpolation::partiallyCoveredSlaveFaces() const +{ + if (!partiallyCoveredSlaveAddrPtr_) + { + calcAddressing(); + } + + return *partiallyCoveredSlaveAddrPtr_; +} + + +template +const scalarField& +GGIInterpolation::slaveFaceUncoveredFractions() const +{ + if (!slaveFaceUncoveredFractionsPtr_) + { + calcAddressing(); + } + + return *slaveFaceUncoveredFractionsPtr_; +} + + template bool GGIInterpolation::movePoints ( diff --git a/src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H b/src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H index 8a5e9e598..282f183de 100644 --- a/src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H +++ b/src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H @@ -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& bridgeField, Field& 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& 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 diff --git a/src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C b/src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C index 48845d702..8b8009735 100644 --- a/src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C +++ b/src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C @@ -68,6 +68,16 @@ GGIInterpolation::featureCosTol_ ); +template +const Foam::debug::tolerancesSwitch +GGIInterpolation::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 @@ -80,7 +90,11 @@ void GGIInterpolation::calcAddressing() const || slaveAddrPtr_ || slaveWeightsPtr_ || uncoveredMasterAddrPtr_ + || partiallyCoveredMasterAddrPtr_ + || masterFaceUncoveredFractionsPtr_ || uncoveredSlaveAddrPtr_ + || partiallyCoveredSlaveAddrPtr_ + || slaveFaceUncoveredFractionsPtr_ ) { FatalErrorIn @@ -601,6 +615,29 @@ void GGIInterpolation::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::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::findNonOverlappingFaces return tpatchFaceNonOverlapAddr; } + +template +void GGIInterpolation::calcPartiallyCoveredFaces +( + const scalarListList& patchWeights, + const scalar& nonOverlapFaceTol, + const bool isMaster +) const +{ + // Sanity checks first + if + ( + isMaster + && (partiallyCoveredMasterAddrPtr_ || masterFaceUncoveredFractionsPtr_) + ) + { + FatalErrorIn + ( + "void GGIInterpolation::" + "calcPartiallyCoveredFaces() const" + ) << "Already calculated master partially covered faces" + << abort(FatalError); + } + + if + ( + !isMaster + && (partiallyCoveredSlaveAddrPtr_ || slaveFaceUncoveredFractionsPtr_) + ) + { + FatalErrorIn + ( + "void GGIInterpolation::" + "calcPartiallyCoveredFaces() const" + ) << "Already calculated slave partially covered faces" + << abort(FatalError); + } + + // Temporary storage + DynamicList patchFacePartialOverlap(patchWeights.size()); + DynamicList 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 void GGIInterpolation:: calcMasterPointAddressing() const diff --git a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C index c9f6bc334..a1d6998e6 100644 --- a/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C +++ b/src/foam/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C @@ -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 );