tried to numerically stabilize geometric function

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 );

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> )

tests/geometry/Functions.cpp

+//======================================================================================================================
+//
+//  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;
+}