diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C index d2a77a26f..3d1b58923 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C @@ -116,11 +116,11 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const forAll (owner, faceI) { - label own = owner[faceI]; - label nei = neighbour[faceI]; + const label own = owner[faceI]; + const label nei = neighbour[faceI]; - vector d = C[nei] - C[own]; - symmTensor wdd = (1.0/magSqr(d))*sqr(d); + const vector d = C[nei] - C[own]; + const symmTensor wdd = (1.0/magSqr(d))*sqr(d); dd[own] += wdd; dd[nei] += wdd; @@ -134,21 +134,10 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010 const vectorField pd = p.delta(); - if (p.coupled()) + forAll (pd, pFaceI) { - forAll (pd, pFaceI) - { - const vector& d = pd[pFaceI]; - dd[fc[pFaceI]] += (1.0/magSqr(d))*sqr(d); - } - } - else - { - forAll (pd, pFaceI) - { - const vector& d = pd[pFaceI]; - dd[fc[pFaceI]] += (1.0/magSqr(d))*sqr(d); - } + const vector& d = pd[pFaceI]; + dd[fc[pFaceI]] += (1.0/magSqr(d))*sqr(d); } } @@ -186,16 +175,20 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const // Least squares vectors on internal faces forAll (owner, faceI) { - label own = owner[faceI]; - label nei = neighbour[faceI]; + const label own = owner[faceI]; + const label nei = neighbour[faceI]; - vector d = C[nei] - C[own]; - scalar magSfByMagSqrd = 1.0/magSqr(d); + const vector d = C[nei] - C[own]; + const scalar magSfByMagSqrd = 1.0/magSqr(d); lsPIn[faceI] = magSfByMagSqrd*(invDd[own] & d); lsNIn[faceI] = -magSfByMagSqrd*(invDd[nei] & d); } + // Get inverse dd boundary field to grab the neighbouring side eventually + const volSymmTensorField::GeometricBoundaryField& volInvDdb = + volInvDd.boundaryField(); + // Least squares vectors on boundary faces forAll (lsP.boundaryField(), patchI) { @@ -204,13 +197,14 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const const fvPatch& p = mesh().boundary()[patchI]; const unallocLabelList& fc = p.faceCells(); + const fvPatchSymmTensorField& volInvDdp = volInvDdb[patchI]; + // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010 const vectorField pd = p.delta(); - if (p.coupled()) + if (volInvDdp.coupled()) { - const symmTensorField invDdNei = - volInvDd.boundaryField()[patchI].patchNeighbourField(); + const symmTensorField invDdNei = volInvDdp.patchNeighbourField(); forAll (pd, pFaceI) { @@ -347,8 +341,8 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const // Internal faces forAll (owner, faceI) { - label own = owner[faceI]; - label nei = neighbour[faceI]; + const label own = owner[faceI]; + const label nei = neighbour[faceI]; if (uggIn[own] > SMALL) { @@ -363,6 +357,10 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const } } + // Get boundary field for cell volumes to fetch the other side eventually + const volScalarField::GeometricBoundaryField& cellVb = + cellV.boundaryField(); + // Boundary faces forAll (lsP.boundaryField(), patchI) { @@ -372,16 +370,16 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const const unallocLabelList& fc = p.faceCells(); const fvsPatchScalarField& pw = w.boundaryField()[patchI]; + const fvPatchScalarField& cellVp = cellVb[patchI]; const vectorField& pSf = Sf.boundaryField()[patchI]; // Get indicator for local side const fvPatchScalarField& ugg = useGaussGrad.boundaryField()[patchI]; const scalarField pUgg = ugg.patchInternalField(); - if (p.coupled()) + if (cellVp.coupled()) { - const scalarField cellVInNei = - cellV.boundaryField()[patchI].patchNeighbourField(); + const scalarField cellVInNei = cellVp.patchNeighbourField(); // Get indicator for neighbour side const scalarField nUgg = ugg.patchNeighbourField(); @@ -417,40 +415,6 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const } } - // Manual check of least squares vectors - - vectorField sumLsq(mesh().nCells(), vector::zero); - vectorField sumSf(mesh().nCells(), vector::zero); - - const vectorField& sfIn = mesh().Sf().internalField(); - - // Least squares vectors on internal faces - forAll (owner, faceI) - { - sumLsq[owner[faceI]] += lsPIn[faceI]; - sumLsq[neighbour[faceI]] += lsNIn[faceI]; - - sumSf[owner[faceI]] += sfIn[faceI]; - sumSf[neighbour[faceI]] -= sfIn[faceI]; - } - - // Least squares vectors on boundary faces - forAll (lsP.boundaryField(), patchI) - { - const vectorField& patchLsP = lsP.boundaryField()[patchI]; - const vectorField& patchSf = mesh().Sf().boundaryField()[patchI]; - - const unallocLabelList& fc = mesh().boundary()[patchI].faceCells(); - - forAll (fc, pFaceI) - { - //sumLsq[fc[pFaceI]] += 0.5*patchLsP[pFaceI]; // works!!! - sumLsq[fc[pFaceI]] += patchLsP[pFaceI]; - - sumSf[fc[pFaceI]] += patchSf[pFaceI]; - } - } - if (debug) { Info<< "leastSquaresVectors::makeLeastSquaresVectors() :" diff --git a/src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H b/src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H index 29e941cef..1d2981d5e 100644 --- a/src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H +++ b/src/foam/interpolations/GGIInterpolation/GGIInterpolationTemplate.H @@ -242,7 +242,7 @@ class GGIInterpolation mutable scalarField* slaveFaceUncoveredFractionsPtr_; - // Private static data + // Private static data //- Facet area error tolerance static const debug::tolerancesSwitch areaErrorTol_; @@ -253,38 +253,35 @@ 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. + // + // Tests using GGI patches of up to ~100K facets have shown that + // the following values gave the best compromise between the + // "octree creation time" vs "octree search time": + // + // octreeSearchMinNLevel_ : 3 + // octreeSearchMaxLeafRatio_ : 3 + // octreeSearchMaxShapeRatio_ : 1 + // + // For GGI patches larger than ~100K facets, your mileage may vary. + // So these 3 control parameters are adjustable using the following + // global optimization switches: + // + // GGIOctreeSearchMinNLevel + // GGIOctreeSearchMaxLeafRatio + // GGIOctreeSearchMaxShapeRatio + // - //- The next 3 attributes are parameters controlling the creation - // of an octree search engine for the GGI facets neighbourhood - // determination. - // - // Tests using GGI patches of up to ~100K facets have shown that - // the following values gave the best compromise between the - // "octree creation time" vs "octree search time": - // - // octreeSearchMinNLevel_ : 3 - // octreeSearchMaxLeafRatio_ : 3 - // octreeSearchMaxShapeRatio_ : 1 - // - // For GGI patches larger than ~100K facets, your mileage may vary. - // So these 3 control parameters are adjustable using the following - // global optimization switches: - // - // GGIOctreeSearchMinNLevel - // GGIOctreeSearchMaxLeafRatio - // GGIOctreeSearchMaxShapeRatio - // + //- Octree search: minNlevel parameter for octree constructor + static const debug::optimisationSwitch octreeSearchMinNLevel_; - //- Octree search: minNlevel parameter for octree constructor - static const debug::optimisationSwitch octreeSearchMinNLevel_; + //- Octree search: maxLeafRatio parameter for octree constructor + static const debug::optimisationSwitch octreeSearchMaxLeafRatio_; - //- Octree search: maxLeafRatio parameter for octree constructor - static const debug::optimisationSwitch octreeSearchMaxLeafRatio_; - - //- Octree search: maxShapeRatio parameter for octree constructor - static const debug::optimisationSwitch octreeSearchMaxShapeRatio_; + //- Octree search: maxShapeRatio parameter for octree constructor + static const debug::optimisationSwitch octreeSearchMaxShapeRatio_; // Private Member Functions @@ -307,6 +304,7 @@ class GGIInterpolation //- Calculate point weights void calcSlavePointWeights() const; + // Helper functions for parallel search //- Is parallel cutting active? @@ -582,6 +580,14 @@ public: ~GGIInterpolation(); + // Public static data + + //- Tolerance for uncovered face areas. Public since the + // MixingPlaneInterpolation needs to be able to change it. Seems more + // appropriate than having friends. + static debug::tolerancesSwitch uncoveredFaceAreaTol_; + + // Member Functions // Access diff --git a/src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C b/src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C index 9d9d41327..f3f80c57b 100644 --- a/src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C +++ b/src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C @@ -69,7 +69,7 @@ GGIInterpolation::featureCosTol_ template -const Foam::debug::tolerancesSwitch +Foam::debug::tolerancesSwitch GGIInterpolation::uncoveredFaceAreaTol_ ( "GGIUncoveredFaceAreaTol", diff --git a/src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationAddressing.C b/src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationAddressing.C index d6908cfc7..da9fbe3ce 100644 --- a/src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationAddressing.C +++ b/src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationAddressing.C @@ -125,6 +125,30 @@ void MixingPlaneInterpolation::calcAddressing() const // introduced by the possibly different // master and slave patch mesh resolution. + // Note: The weights in the GGI interpolation must be scaled when using + // interpolation from/to ribbon patches. Now, the GGI interpolation can also + // take care of partially overlapping faces where we must not scale the + // weights (e.g. bridgeOverlap option for partial symmetry plane + // treatment in uncovered faces). In order to avoid having partially + // overlapping faces in the mixing plane, we will manually set the + // corresponding tolerance here and reset it after we calculate the + // addressing. VV, 31/Aug/2018. + // Save old tolerance + const scalar oldTolerance = + GGIInterpolation + < + standAlonePatch, + standAlonePatch + >::uncoveredFaceAreaTol_(); + + // Set tolerance to zero such that no faces qualify as partially + // overlapping faces and all the weights are properly scaled + GGIInterpolation + < + standAlonePatch, + standAlonePatch + >::uncoveredFaceAreaTol_ = 0.0; + GGIInterpolation masterCircumAvgPatchToPatch ( @@ -202,6 +226,16 @@ void MixingPlaneInterpolation::calcAddressing() const slaveProfileToPatchWeightsPtr_ = new scalarListList(slaveCircumAvgPatchToPatch.slaveWeights()); + + + // Now that we have finished calculating the addressing and the weights, + // let's reset the old tolerances in case some other GGIInterpolation with + // stand alone patches uses the partially covered face machinery. + GGIInterpolation + < + standAlonePatch, + standAlonePatch + >::uncoveredFaceAreaTol_ = oldTolerance; } diff --git a/src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationTemplate.H b/src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationTemplate.H index 3380f31c6..128962d26 100644 --- a/src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationTemplate.H +++ b/src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationTemplate.H @@ -53,11 +53,11 @@ Contributor Giving credit where credit is due: 1: Hakan Nilsson from the Chalmers University of Technology came up with - the initial idea of using a combination of 2 GGI interfaces sharing a - common single-face 360 degree ribbons patch in order to compute the - circumferential average of fields. This implementation of the - mixingPlane interpolation algorithm is an exploration of this simple - but rather powerful idea. + the initial idea of using a combination of 2 GGI interfaces sharing a + common single-face 360 degree ribbons patch in order to compute the + circumferential average of fields. This implementation of the + mixingPlane interpolation algorithm is an exploration of this simple + but rather powerful idea. 2: Maryse Page from Hydro-Quebec provided many test cases and many simulation runs for testing this interpolation algorithm. Testing is @@ -240,7 +240,7 @@ class MixingPlaneInterpolation //- Interpolation profile // This list of points defines the profile generating the 'n' // circular bands we are going to use for - // the circumferential averaging algorithm. For for 'n' + // the circumferential averaging algorithm. For 'n' // circular bands, we need 'n + 1' points. mutable pointField interpolationProfile_; diff --git a/src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneProfile.C b/src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneProfile.C index 18cc00822..ebc342db5 100644 --- a/src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneProfile.C +++ b/src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneProfile.C @@ -313,13 +313,13 @@ MixingPlaneInterpolation::calcProfile() const if (debug > 1) { // Write histograms - forAllIter (profileHistogram, masterHistogram, zHi) + forAllConstIter (profileHistogram, masterHistogram, zHi) { Info<< "master histo (z, n): (" << zHi->first << " " << static_cast(zHi->second.size()) << ")" << endl; } - forAllIter (profileHistogram, slaveHistogram, zHi) + forAllConstIter (profileHistogram, slaveHistogram, zHi) { Info<< "slave histo (z, n): (" << zHi->first << " " << static_cast(zHi->second.size()) << ")" << endl; diff --git a/src/foam/meshes/polyMesh/polyMesh.C b/src/foam/meshes/polyMesh/polyMesh.C index 6ef9def6d..21f1d3f5f 100644 --- a/src/foam/meshes/polyMesh/polyMesh.C +++ b/src/foam/meshes/polyMesh/polyMesh.C @@ -1232,10 +1232,11 @@ Foam::tmp Foam::polyMesh::movePoints geometricD_ = Vector