Merge pull request #255 from MFraters/refactor_and_fix_spherical_slabs

Refactor and fix spherical slabs
GeodynamicWorldBuilder · Jun 12, 2021 · b3eabfb · b3eabfb
2 parents 7fa5556 + 35cd398
commit b3eabfb
Show file tree

Hide file tree

Showing 10 changed files with 447 additions and 189 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -25,6 +25,7 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm
 - The Vizualizer now uses compressed output by default. This decreases the file size and increases perforamnce. \[Menno Fraters; 2021-05-22; [#239](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/239)\]
 - The Vizualizer buffers output before it writes it to a file. This increases performance. \[Menno Fraters; 2021-05-22; [#239](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/239)\]
 - In the fault temperature model linear, the entiry top temperature is changed to center temperature and the entry bottom temperature is changed to side temperature, since this represents the actual sides more accurately. \[Menno Fraters; 2021-07-09; [#260](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/260)\]
+- A large overhoal of the distance_point_from_curved_planes function improving the accuracy in spherical coordinates. This may slighty change the results, but it should be an improvement both in accuracy and performance. Also makes available a coordinate system independent distrance function and some fast trigonometric functions. \[Menno Fraters; 2021-06-11; [#255](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/255)\]
 
 ### Fixed
 - Fixed namespaces, adding `WorldBuilder::` where needed. \[Timo Heister; 2020-08-10; [#205](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/205)\] 

diff --git a/include/world_builder/point.h b/include/world_builder/point.h
@@ -25,6 +25,7 @@
 #include <iostream>
 
 #include <world_builder/coordinate_system.h>
+#include <world_builder/assert.h>
 
 namespace WorldBuilder
 {
@@ -133,14 +134,40 @@ namespace WorldBuilder
       /**
        * access index (const)
        */
-      const double &operator[](const unsigned int index) const;
+      inline
+      const double &operator[](const unsigned int index) const
+      {
+        WBAssert(index <= dim, "Can't ask for element " << index << " in an point with dimension " << dim << ".");
+        return point[index];
+      }
 
 
       /**
        * access index
        */
-      double &operator[](const unsigned int index);
+      inline double &operator[](const unsigned int index)
+      {
+        WBAssert(index <= dim, "Can't ask for element " << index << " in an point with dimension " << dim << ".");
+        return point[index];
+      }
+
+      /**
+       * Computes the distance between this and a given point.
+       * In spherical coordinates it returns the central angle in radians.
+       */
+      double
+      distance(const Point<dim> &two) const;
 
+      /**
+       * Computes the cheapest relative distance between this and a given point.
+       * The return value itself is only guartenteed to have the property that a
+       * larger value is further away.
+       * In the current implementation that means for the cartasian case the squared
+       * value is returned and for the spherical value the result of the havearsine
+       * function without asin and sqrt is returned.
+       */
+      double
+      cheap_relative_distance(const Point<dim> &two) const;
 
       /**
        * return the internal array which stores the point data.
@@ -188,6 +215,73 @@ namespace WorldBuilder
 
   };
 
+  /**
+   * This namespace contains some faster but less accurate version of the
+   * trigonomic functions and a faster version of the fmod function.
+   */
+  namespace FT
+  {
+    constexpr double const_pi = 3.141592653589793238462643383279502884;
+
+    /**
+     * Fast version of the fmod function.
+     */
+    inline double fmod(const double x, const double y)
+    {
+      const double x_div_y = x/y;
+      return (x_div_y-(int)x_div_y)*y;
+    }
+
+    /**
+     * Fast sin function, accurate for values between 0 and pi. The implemenation is
+     * based on discussion at https://stackoverflow.com/a/6104692.
+     *
+     * The accuracy seem good enough for most purposes. The unit test tests in steps
+     * of 0.01 from -4 pi to 4 pi and compares against the std sin function and the difference
+     * is always smaller than 1.2e-5. If the test is run with intervals of 0.001 then there
+     * are 12 entries which are (very slightly) above that (<3e-8) at angles of about
+     * -174, -6, 6  and 174.
+     *
+     */
+    inline double fast_sin_d(const double angle)
+    {
+      constexpr double A = 4.0/(const_pi *const_pi);
+      constexpr double oneminPmin = 1.-0.1952403377008734-0.01915214119105392;
+
+      const double y = A* angle * ( const_pi - angle );
+      return y*( oneminPmin + y*( 0.1952403377008734 + y * 0.01915214119105392 ) ) ;
+    }
+
+    /**
+     * Fast but less accurate sin function for any angle.
+     * Implemented by calling fast_sin_d with a mirrored x if needed to
+     * forfill the constrained of fast_sin_d to only have values between
+     * zero and pi.
+     */
+    inline double sin(const double raw_angle)
+    {
+      const double angle = (raw_angle > -const_pi && raw_angle < const_pi)
+                           ?
+                           raw_angle
+                           :
+                           FT::fmod(raw_angle + std::copysign(const_pi,raw_angle), const_pi * 2.0) - std::copysign(const_pi,raw_angle);
+
+      if (angle >= 0)
+        return fast_sin_d(angle);
+      else
+        return -fast_sin_d(-angle);
+    }
+
+    /**
+     * Fast but less accurate cos function for any angle.
+     */
+    inline double cos(const double angle)
+    {
+      return FT::sin((const_pi*0.5)-angle);
+    }
+  }
+
+
   template<int dim>
   Point<dim> operator*(const double scalar, const Point<dim> &point);
 

diff --git a/include/world_builder/utilities.h b/include/world_builder/utilities.h
@@ -232,7 +232,7 @@ namespace WorldBuilder
         /**
          * Evaluate at point @p x.
          */
-        double operator() (double x) const;
+        double operator() (const double x) const;
 
       private:
         /**

diff --git a/source/point.cc b/source/point.cc
@@ -205,30 +205,6 @@ namespace WorldBuilder
   }
 
 
-  /**
-   * access index
-   */
-  template<int dim>
-  const double &
-  Point<dim>::operator[](const unsigned int index) const
-  {
-    WBAssertThrow(index <= dim, "Can't ask for element " << index << " in an point with dimension " << dim << ".");
-    return point[index];
-  }
-
-
-  /**
-   * access index
-   */
-  template<int dim>
-  double &
-  Point<dim>::operator[](const unsigned int index)
-  {
-    WBAssertThrow(index <= dim, "Can't ask for element " << index << " in an point with dimension " << dim << ".");
-    return point[index];
-  }
-
-
   template<int dim>
   const std::array<double,dim> &
   Point<dim>::get_array() const
@@ -268,6 +244,52 @@ namespace WorldBuilder
   }
 
 
+  template<int dim>
+  double
+  Point<dim>::distance(const Point<dim> &two) const
+  {
+    if (this->coordinate_system == spherical)
+      {
+        // spherical
+        const double d_longitude = two[0] - this->point[0];
+        const double d_lattitude = two[1] - this->point[1];
+        const double sin_d_lat = std::sin(d_lattitude * 0.5);
+        const double sin_d_long = std::sin(d_longitude * 0.5);
+        return 2.0 * asin(sqrt((sin_d_lat * sin_d_lat) + (sin_d_long*sin_d_long) * std::cos(this->point[1]) * std::cos(two[1])));
+      }
+    else
+      {
+        // cartesian
+        const double x_distance_to_reference_point = point[0]-two[0];
+        const double y_distance_to_reference_point = point[1]-two[1];
+        return sqrt((x_distance_to_reference_point*x_distance_to_reference_point) + (y_distance_to_reference_point*y_distance_to_reference_point));
+      }
+  }
+
+
+
+  template<int dim>
+  double
+  Point<dim>::cheap_relative_distance(const Point<dim> &two) const
+  {
+    if (this->coordinate_system == spherical)
+      {
+        // spherical
+        const double d_longitude = two[0] - this->point[0];
+        const double d_lattitude = two[1] - this->point[1];
+        const double sin_d_lat = FT::sin(d_lattitude * 0.5);
+        const double sin_d_long = FT::sin(d_longitude * 0.5);
+        return (sin_d_lat * sin_d_lat) + (sin_d_long*sin_d_long) * FT::cos(this->point[1]) * FT::cos(two[1]);
+      }
+    else
+      {
+        // cartesian
+        const double x_distance_to_reference_point = point[0]-two[0];
+        const double y_distance_to_reference_point = point[1]-two[1];
+        return (x_distance_to_reference_point*x_distance_to_reference_point) + (y_distance_to_reference_point*y_distance_to_reference_point);
+      }
+  }
+
   /**
    * Multiplies a point with a scalar.
    */