exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Added stabilisation in division by zero; optimising compiler

Browse source
This commit is contained in:
Hrvoje Jasak 2011-09-22 11:30:24 +01:00
parent 63f0794020
commit 09b8e88887
1 changed files with 38 additions and 8 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

46
src/OpenFOAM/algorithms/octree/octree/treeBoundBox.C
View file

@ -318,8 +318,21 @@ bool Foam::treeBoundBox::intersects
) const
{
const direction endBits = posBits(end);
pt = start;
// Note
// Optimising compilers (Intel 12) reorganises the code in
// a way which executes all divisions before doing the if-check.
// This causes a floating point exception of division by zero.
// Furthermore, division by VSMALL (first number larger than zero) may
// in some cases cause floating point exception overflow, and stabilisation
// is required. Added by HJ, 21/Sep/2011
scalar s;
static const scalar smallNumber = SMALL;
while (true)
{
direction ptBits = posBits(pt);
@ -328,6 +341,7 @@ bool Foam::treeBoundBox::intersects
{
// pt inside bb
ptOnFaces = faceBits(pt);
return true;
}
@ -335,6 +349,7 @@ bool Foam::treeBoundBox::intersects
{
// pt and end in same block outside of bb
ptOnFaces = faceBits(pt);
return false;
}
@ -343,7 +358,9 @@ bool Foam::treeBoundBox::intersects
// Intersect with plane V=min, n=-1,0,0
if (Foam::mag(overallVec.x()) > VSMALL)
{
scalar s = (min().x() - overallStart.x())/overallVec.x();
s = (min().x() - overallStart.x())/
stabilise(overallVec.x(), smallNumber);
pt.x() = min().x();
pt.y() = overallStart.y() + overallVec.y()*s;
pt.z() = overallStart.z() + overallVec.z()*s;
@ -360,7 +377,9 @@ bool Foam::treeBoundBox::intersects
// Intersect with plane V=max, n=1,0,0
if (Foam::mag(overallVec.x()) > VSMALL)
{
scalar s = (max().x() - overallStart.x())/overallVec.x();
s = (max().x() - overallStart.x())/
stabilise(overallVec.x(), smallNumber);
pt.x() = max().x();
pt.y() = overallStart.y() + overallVec.y()*s;
pt.z() = overallStart.z() + overallVec.z()*s;
@ -375,7 +394,9 @@ bool Foam::treeBoundBox::intersects
// Intersect with plane V=min, n=0,-1,0
if (Foam::mag(overallVec.y()) > VSMALL)
{
scalar s = (min().y() - overallStart.y())/overallVec.y();
s = (min().y() - overallStart.y())/
stabilise(overallVec.y(), smallNumber);
pt.x() = overallStart.x() + overallVec.x()*s;
pt.y() = min().y();
pt.z() = overallStart.z() + overallVec.z()*s;
@ -390,7 +411,9 @@ bool Foam::treeBoundBox::intersects
// Intersect with plane V=max, n=0,1,0
if (Foam::mag(overallVec.y()) > VSMALL)
{
scalar s = (max().y() - overallStart.y())/overallVec.y();
s = (max().y() - overallStart.y())/
stabilise(overallVec.y(), smallNumber);
pt.x() = overallStart.x() + overallVec.x()*s;
pt.y() = max().y();
pt.z() = overallStart.z() + overallVec.z()*s;
@ -405,7 +428,8 @@ bool Foam::treeBoundBox::intersects
// Intersect with plane V=min, n=0,0,-1
if (Foam::mag(overallVec.z()) > VSMALL)
{
scalar s = (min().z() - overallStart.z())/overallVec.z();
s = (min().z() - overallStart.z())/
stabilise(overallVec.z(), smallNumber);
pt.x() = overallStart.x() + overallVec.x()*s;
pt.y() = overallStart.y() + overallVec.y()*s;
pt.z() = min().z();
@ -420,7 +444,9 @@ bool Foam::treeBoundBox::intersects
// Intersect with plane V=max, n=0,0,1
if (Foam::mag(overallVec.z()) > VSMALL)
{
scalar s = (max().z() - overallStart.z())/overallVec.z();
s = (max().z() - overallStart.z())/
stabilise(overallVec.z(), smallNumber);
pt.x() = overallStart.x() + overallVec.x()*s;
pt.y() = overallStart.y() + overallVec.y()*s;
pt.z() = max().z();
@ -442,7 +468,8 @@ bool Foam::treeBoundBox::intersects
) const
{
direction ptBits;
return intersects(start, end-start, start, end, pt, ptBits);
return intersects(start, end - start, start, end, pt, ptBits);
}
@ -458,7 +485,7 @@ bool Foam::treeBoundBox::contains(const vector& dir, const point& pt) const
//
// Compare all components against min and max of bb
//
for (direction cmpt=0; cmpt<3; cmpt++)
for (direction cmpt = 0; cmpt < 3; cmpt++)
{
if (pt[cmpt] < min()[cmpt])
{
@ -496,6 +523,7 @@ bool Foam::treeBoundBox::contains(const vector& dir, const point& pt) const
Foam::direction Foam::treeBoundBox::faceBits(const point& pt) const
{
direction faceBits = 0;
if (pt.x() == min().x())
{
faceBits |= LEFTBIT;
@ -522,6 +550,7 @@ Foam::direction Foam::treeBoundBox::faceBits(const point& pt) const
{
faceBits |= FRONTBIT;
}
return faceBits;
}
@ -557,6 +586,7 @@ Foam::direction Foam::treeBoundBox::posBits(const point& pt) const
{
posBits |= FRONTBIT;
}
return posBits;
}