From 559e9a8f681d26a3a4c529cf2ba47926972215c9 Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Mon, 10 Jul 2017 10:35:57 +0200
Subject: [PATCH] tried to numerically stabilize geometric function
---
src/geometry/GeometricalFunctions.cpp | 26 +++++++-------
tests/geometry/CMakeLists.txt | 2 ++
tests/geometry/Functions.cpp | 52 +++++++++++++++++++++++++++
3 files changed, 67 insertions(+), 13 deletions(-)
create mode 100644 tests/geometry/Functions.cpp
diff --git a/src/geometry/GeometricalFunctions.cpp b/src/geometry/GeometricalFunctions.cpp
index 52f3c7cb9..ec8e02b82 100644
--- a/src/geometry/GeometricalFunctions.cpp
+++ b/src/geometry/GeometricalFunctions.cpp
@@ -259,7 +259,7 @@ void getClosestLineSegmentPoints( const Vector3<real_t>& a1, const Vector3<real_
const Vector3<real_t> a1b1( b1 - a1 );
const real_t da1 ( a1a2 * a1b1 );
const real_t db1 ( b1b2 * a1b1 );
- if( da1 <= real_t(0) && db1 >= real_t(0) ) {
+ if( da1 <= +Limits<real_t>::fpuAccuracy() && db1 >= -Limits<real_t>::fpuAccuracy() ) {
cp1 = a1;
cp2 = b1;
return;
@@ -268,7 +268,7 @@ void getClosestLineSegmentPoints( const Vector3<real_t>& a1, const Vector3<real_
const Vector3<real_t> a1b2( b2 - a1 );
const real_t da2 ( a1a2 * a1b2 );
const real_t db2 ( b1b2 * a1b2 );
- if( da2 <= real_t(0) && db2 <= real_t(0) ) {
+ if( da2 <= +Limits<real_t>::fpuAccuracy() && db2 <= +Limits<real_t>::fpuAccuracy() ) {
cp1 = a1;
cp2 = b2;
return;
@@ -277,7 +277,7 @@ void getClosestLineSegmentPoints( const Vector3<real_t>& a1, const Vector3<real_
const Vector3<real_t> a2b1( b1 - a2 );
const real_t da3 ( a1a2 * a2b1 );
const real_t db3 ( b1b2 * a2b1 );
- if( da3 >= real_t(0) && db3 >= real_t(0) ) {
+ if( da3 >= -Limits<real_t>::fpuAccuracy() && db3 >= -Limits<real_t>::fpuAccuracy() ) {
cp1 = a2;
cp2 = b1;
return;
@@ -286,7 +286,7 @@ void getClosestLineSegmentPoints( const Vector3<real_t>& a1, const Vector3<real_
const Vector3<real_t> a2b2( b2 - a2 );
const real_t da4 ( a1a2 * a2b2 );
const real_t db4 ( b1b2 * a2b2 );
- if( da4 >= real_t(0) && db4 <= real_t(0) ) {
+ if( da4 >= -Limits<real_t>::fpuAccuracy() && db4 <= +Limits<real_t>::fpuAccuracy() ) {
cp1 = a2;
cp2 = b2;
return;
@@ -300,20 +300,20 @@ void getClosestLineSegmentPoints( const Vector3<real_t>& a1, const Vector3<real_
Vector3<real_t> n, k;
const real_t la( a1a2 * a1a2 );
- if( da1 >= real_t(0) && da3 <= real_t(0) ) {
+ if( da1 >= -Limits<real_t>::fpuAccuracy() && da3 <= +Limits<real_t>::fpuAccuracy() ) {
k = (da1 / la) * a1a2;
n = a1b1 - k;
- if( b1b2 * n >= real_t(0) ) {
+ if( b1b2 * n >= -Limits<real_t>::fpuAccuracy() ) {
cp1 = a1 + k;
cp2 = b1;
return;
}
}
- if( da2 >= real_t(0) && da4 <= real_t(0) ) {
+ if( da2 >= -Limits<real_t>::fpuAccuracy() && da4 <= +Limits<real_t>::fpuAccuracy() ) {
k = (da2 / la) * a1a2;
n = a1b2 - k;
- if( b1b2 * n <= real_t(0) ) {
+ if( b1b2 * n <= +Limits<real_t>::fpuAccuracy() ) {
cp1 = a1 + k;
cp2 = b2;
return;
@@ -321,20 +321,20 @@ void getClosestLineSegmentPoints( const Vector3<real_t>& a1, const Vector3<real_
}
const real_t lb( b1b2 * b1b2 );
- if( db1 <= real_t(0) && db2 >= real_t(0) ) {
+ if( db1 <= +Limits<real_t>::fpuAccuracy() && db2 >= -Limits<real_t>::fpuAccuracy() ) {
k = (-db1 / lb) * b1b2;
n = -a1b1 - k;
- if( a1a2 * n >= real_t(0) ) {
+ if( a1a2 * n >= -Limits<real_t>::fpuAccuracy() ) {
cp1 = a1;
cp2= b1 + k;
return;
}
}
- if( db3 <= real_t(0) && db4 >= real_t(0) ) {
+ if( db3 <= +Limits<real_t>::fpuAccuracy() && db4 >= -Limits<real_t>::fpuAccuracy() ) {
k = (-db3 / lb) * b1b2;
n = -a2b1 - k;
- if( a1a2 * n <= real_t(0) ) {
+ if( a1a2 * n <= +Limits<real_t>::fpuAccuracy() ) {
cp1 = a2;
cp2 = b1 + k;
return;
@@ -347,7 +347,7 @@ void getClosestLineSegmentPoints( const Vector3<real_t>& a1, const Vector3<real_
const real_t scale( a1a2 * b1b2 );
real_t div( la * lb - math::sq(scale) );
- WALBERLA_ASSERT_GREATER( div, real_t(0), "Invalid division" );
+ WALBERLA_ASSERT_GREATER( div, real_t(0), std::setprecision(16) << "Invalid division\n" << a1 << "\n" << a2 << "\n" << b1 << "\n" << b2 );
div = real_t(1) / div;
const real_t s( ( lb * da1 - scale * db1 ) * div );
diff --git a/tests/geometry/CMakeLists.txt b/tests/geometry/CMakeLists.txt
index 2bedbd39e..789e44395 100644
--- a/tests/geometry/CMakeLists.txt
+++ b/tests/geometry/CMakeLists.txt
@@ -8,6 +8,8 @@
waLBerla_compile_test( FILES CellAABBTest.cpp )
waLBerla_execute_test( NAME CellAABBTest COMMAND $<TARGET_FILE:CellAABBTest> )
+waLBerla_compile_test( FILES Functions.cpp )
+waLBerla_execute_test( NAME Functions )
waLBerla_compile_test( FILES VoxelFileTest.cpp )
waLBerla_execute_test( NAME VoxelFileTest COMMAND $<TARGET_FILE:VoxelFileTest> )
diff --git a/tests/geometry/Functions.cpp b/tests/geometry/Functions.cpp
new file mode 100644
index 000000000..f4341fae6
--- /dev/null
+++ b/tests/geometry/Functions.cpp
@@ -0,0 +1,52 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file Functions.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include <core/DataTypes.h>
+#include <core/debug/TestSubsystem.h>
+#include <core/logging/Logging.h>
+#include <core/mpi/MPIManager.h>
+#include <geometry/GeometricalFunctions.h>
+
+using namespace walberla;
+using namespace walberla::geometry;
+
+typedef math::Vector3<real_t> Vec3;
+
+int main( int argc, char** argv )
+{
+ debug::enterTestMode();
+
+ MPIManager::instance()->initializeMPI( &argc, &argv );
+
+ Vec3 a,b;
+ geometry::getClosestLineSegmentPoints(Vec3(real_c(7.41657),real_c( 9.82461),real_c(17.28197)),
+ Vec3(real_c(7.21657),real_c( 9.82461),real_c(17.28197)),
+ Vec3(real_c(7.36474),real_c(10.20490),real_c(17.45072)),
+ Vec3(real_c(7.16474),real_c(10.20490),real_c(17.45072)), a, b);
+ WALBERLA_LOG_DEVEL(a << "\n" << b);
+ // might fail due to numerically instable algorithm
+ geometry::getClosestLineSegmentPoints(Vec3(real_c(7.416573195643714),real_c(9.82461330720222),real_c(17.28197302209214)),
+ Vec3(real_c(7.216573195643715),real_c(9.82461330720222),real_c(17.28197302209214)),
+ Vec3(real_c(7.364742740420613),real_c(10.2049088375177),real_c(17.45072353608904)),
+ Vec3(real_c(7.164742740420614),real_c(10.2049088375177),real_c(17.45072353608904)), a, b);
+ WALBERLA_LOG_DEVEL(a << "\n" << b);
+
+ return EXIT_SUCCESS;
+}
--
GitLab