From 9ccc2d8fd5df41c7c8e2e7ac9cd03ae08b69a283 Mon Sep 17 00:00:00 2001 From: Kutalmis Bercin Date: Tue, 12 Jul 2022 16:00:34 +0100 Subject: [PATCH] ENH: edgeInterpolation: avoid division-by-zero errors in skew corrections --- .../edgeInterpolation/edgeInterpolation.C | 25 ++++++++++++++----- 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C index 370e95ff6d..aa44dd3fae 100644 --- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C +++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2016-2017 Wikki Ltd - Copyright (C) 2020 OpenCFD Ltd. + Copyright (C) 2020-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -581,8 +581,15 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const const vector& S = points[edges[edgeI].start()]; vector e = edges[edgeI].vec(points); - scalar alpha = - -(((N - P)^(S - P))&((N - P)^e ))/(((N - P)^e)&((N - P)^e)); + const scalar beta = magSqr((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); @@ -614,9 +621,15 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const const vector& S = points[patchEdges[edgeI].start()]; vector e = patchEdges[edgeI].vec(points); - scalar alpha = - - (((N - P)^(S - P))&((N - P)^e)) - /(((N - P)^e)&((N - P)^e)); + const scalar beta = magSqr((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);