ENH: edgeInterpolation: avoid division-by-zero errors in skew corrections

This commit is contained in:
Kutalmis Bercin
2022-07-12 16:00:34 +01:00
parent c2cae92fc5
commit 9ccc2d8fd5

View File

@ -6,7 +6,7 @@
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd Copyright (C) 2016-2017 Wikki Ltd
Copyright (C) 2020 OpenCFD Ltd. Copyright (C) 2020-2022 OpenCFD Ltd.
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -581,8 +581,15 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
const vector& S = points[edges[edgeI].start()]; const vector& S = points[edges[edgeI].start()];
vector e = edges[edgeI].vec(points); vector e = edges[edgeI].vec(points);
scalar alpha = const scalar beta = magSqr((N - P)^e);
-(((N - P)^(S - P))&((N - P)^e ))/(((N - P)^e)&((N - P)^e));
if (beta < ROOTVSMALL)
{
// Too small - skew correction remains zero
continue;
}
const scalar alpha = -(((N - P)^(S - P))&((N - P)^e))/beta;
vector E(S + alpha*e); vector E(S + alpha*e);
@ -614,9 +621,15 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
const vector& S = points[patchEdges[edgeI].start()]; const vector& S = points[patchEdges[edgeI].start()];
vector e = patchEdges[edgeI].vec(points); vector e = patchEdges[edgeI].vec(points);
scalar alpha = const scalar beta = magSqr((N - P)^e);
- (((N - P)^(S - P))&((N - P)^e))
/(((N - P)^e)&((N - P)^e)); if (beta < ROOTVSMALL)
{
// Too small - skew correction remains zero
continue;
}
const scalar alpha = -(((N - P)^(S - P))&((N - P)^e))/beta;
vector E(S + alpha*e); vector E(S + alpha*e);