From d5b7200295fba538046ff07f94ca653cafb5c26a Mon Sep 17 00:00:00 2001
From: Kutalmis Bercin <kutalmis.bercin@esi-group.com>
Date: Fri, 2 Dec 2022 15:00:15 +0000
Subject: [PATCH] ENH: edgeInterpolation: make skewness-correction optimisation
compile-time optional
---
.../edgeInterpolation/edgeInterpolation.C | 105 +++++++++++-------
.../edgeInterpolation/edgeInterpolation.H | 4 +
2 files changed, 66 insertions(+), 43 deletions(-)
diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C
index aa44dd3faea..ad65518493a 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C
+++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C
@@ -174,6 +174,25 @@ bool Foam::edgeInterpolation::movePoints() const
}
+const Foam::vector& Foam::edgeInterpolation::skewCorr(const label edgeI) const
+{
+ #ifdef FA_SKEW_CORRECTION
+
+ return
+ (
+ skewCorrectionVectors_
+ ? (*skewCorrectionVectors_)[edgeI]
+ : pTraits<vector>::zero
+ );
+
+ #else
+
+ return (*skewCorrectionVectors_)[edgeI];
+
+ #endif
+}
+
+
void Foam::edgeInterpolation::makeLPN() const
{
DebugInFunction
@@ -205,20 +224,18 @@ void Foam::edgeInterpolation::makeLPN() const
scalarField& lPNIn = lPN.primitiveFieldRef();
+ // Calculate skewness correction vectors if necessary
+ (void) skew();
+
forAll(owner, edgeI)
{
- vector curSkewCorrVec(Zero);
-
- if (skew())
- {
- curSkewCorrVec = skewCorrectionVectors()[edgeI];
- }
+ const vector& skewCorrEdge = skewCorr(edgeI);
scalar lPE =
mag
(
edgeCentres[edgeI]
- - curSkewCorrVec
+ - skewCorrEdge
- faceCentres[owner[edgeI]]
);
@@ -227,7 +244,7 @@ void Foam::edgeInterpolation::makeLPN() const
(
faceCentres[neighbour[edgeI]]
- edgeCentres[edgeI]
- + curSkewCorrVec
+ + skewCorrEdge
);
lPNIn[edgeI] = (lPE + lEN);
@@ -283,20 +300,18 @@ void Foam::edgeInterpolation::makeWeights() const
scalarField& weightingFactorsIn = weightingFactors.primitiveFieldRef();
+ // Calculate skewness correction vectors if necessary
+ (void) skew();
+
forAll(owner, edgeI)
{
- vector curSkewCorrVec(Zero);
-
- if (skew())
- {
- curSkewCorrVec = skewCorrectionVectors()[edgeI];
- }
+ const vector& skewCorrEdge = skewCorr(edgeI);
scalar lPE =
mag
(
edgeCentres[edgeI]
- - curSkewCorrVec
+ - skewCorrEdge
- faceCentres[owner[edgeI]]
);
@@ -305,7 +320,7 @@ void Foam::edgeInterpolation::makeWeights() const
(
faceCentres[neighbour[edgeI]]
- edgeCentres[edgeI]
- + curSkewCorrVec
+ + skewCorrEdge
);
weightingFactorsIn[edgeI] =
@@ -370,6 +385,8 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const
const edgeList& edges = mesh().edges();
const pointField& points = mesh().points();
+ // Calculate skewness correction vectors if necessary
+ (void) skew();
forAll(owner, edgeI)
{
@@ -386,19 +403,13 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const
unitDelta.normalise();
- // Calc PN arc length
- vector curSkewCorrVec(Zero);
-
- if (skew())
- {
- curSkewCorrVec = skewCorrectionVectors()[edgeI];
- }
+ const vector& skewCorrEdge = skewCorr(edgeI);
scalar lPE =
mag
(
edgeCentres[edgeI]
- - curSkewCorrVec
+ - skewCorrEdge
- faceCentres[owner[edgeI]]
);
@@ -407,7 +418,7 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const
(
faceCentres[neighbour[edgeI]]
- edgeCentres[edgeI]
- + curSkewCorrVec
+ + skewCorrEdge
);
scalar lPN = lPE + lEN;
@@ -643,46 +654,54 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
}
}
+ #ifdef FA_SKEW_CORRECTION
- scalar skewCoeff = 0.0;
-
- // Calculating PN arc length
- scalarField lPN(owner.size());
+ constexpr scalar maxSkewRatio = 0.1;
+ scalar skewCoeff = 0;
- forAll(owner, edgeI)
+ forAll(own, edgeI)
{
- lPN[edgeI] =
+ const scalar magSkew = mag(skewCorrVecs[edgeI]);
+
+ const scalar lPN =
mag
(
Ce[edgeI]
- - SkewCorrVecs[edgeI]
+ - skewCorrVecs[edgeI]
- C[owner[edgeI]]
)
+ mag
(
C[neighbour[edgeI]]
- Ce[edgeI]
- + SkewCorrVecs[edgeI]
+ + skewCorrVecs[edgeI]
);
- }
- if (lPN.size() > 0)
- {
- skewCoeff = max(mag(SkewCorrVecs.internalField())/mag(lPN));
+ const scalar ratio = magSkew/lPN;
+
+ if (skewCoeff < ratio)
+ {
+ skewCoeff = ratio;
+
+ if (skewCoeff > maxSkewRatio)
+ {
+ break;
+ }
+ }
}
DebugInFunction
<< "skew coefficient = " << skewCoeff << endl;
- if (skewCoeff < 0.1)
+ if (skewCoeff < maxSkewRatio)
{
- skew_ = false;
deleteDemandDrivenData(skewCorrectionVectors_);
}
- else
- {
- skew_ = true;
- }
+
+ #endif
+
+ skew_ = bool(skewCorrectionVectors_);
+
DebugInFunction
<< "Finished constructing skew correction vectors"
diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.H b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.H
index b740b0a07f9..ff84b3c7420 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.H
+++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.H
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd
+ Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -93,6 +94,9 @@ class edgeInterpolation
// Private Member Functions
+ //- Return skewness correction per edge
+ const vector& skewCorr(const label edgeI) const;
+
//- Construct geodesic distance between P and N
void makeLPN() const;
--
GitLab