From 05508d97bd6d1546eafe1cae50c36560a6e6ba9a Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Fri, 12 Jun 2020 10:04:59 +0200
Subject: [PATCH 1/2] [md][embedded] Add optimized version of circle points

---
 dumux/multidomain/embedded/circlepoints.hh | 77 +++++++++++++++-------
 1 file changed, 53 insertions(+), 24 deletions(-)

diff --git a/dumux/multidomain/embedded/circlepoints.hh b/dumux/multidomain/embedded/circlepoints.hh
index 64319d8292..ae9aee1a70 100644
--- a/dumux/multidomain/embedded/circlepoints.hh
+++ b/dumux/multidomain/embedded/circlepoints.hh
@@ -28,38 +28,46 @@
 #include <vector>
 #include <cmath>
 
+#include <dumux/common/math.hh>
 #include <dune/common/exceptions.hh>
 
-namespace Dumux {
-namespace EmbeddedCoupling {
+namespace Dumux::EmbeddedCoupling {
 
 /*!
  * \ingroup EmbeddedCoupling
- * \brief returns a vector of points on a circle
+ * \param points a vector of points to be filled
+ * \param sincos vector with [sin(a0), cos(a0), sin(a1), cos(a1), ...] for each circumferential sample point ai
  * \param center the circle center
  * \param normal the normal to the circle plane
  * \param radius the circle radius
- * \param numPoints the number of points
+ * \note This version allows for a more efficient circle point generator
+ *       if the circumferential positions are fixed and can be reused.
+ *       In this case, sine and cosine of the corresponding angles can be precomputed.
  */
-template<class GlobalPosition>
-std::vector<GlobalPosition> circlePoints(const GlobalPosition& center,
-                                         const GlobalPosition& normal,
-                                         const typename GlobalPosition::value_type radius,
-                                         const std::size_t numPoints = 20)
+template<class GlobalPosition, class Scalar>
+void circlePoints(std::vector<GlobalPosition>& points,
+                  const std::vector<Scalar>& sincos,
+                  const GlobalPosition& center,
+                  const GlobalPosition& normal,
+                  const Scalar radius)
 {
-    using std::abs; using std::sin; using std::cos;
+    assert(sincos.size() % 2 == 0 && "Sample angles have to be pairs of sin/cos so size needs to be even.");
+    static_assert(GlobalPosition::dimension == 3, "Only implemented for world dimension 3");
+
+    using std::abs;
     using ctype = typename GlobalPosition::value_type;
 
     constexpr ctype eps = 1.5e-7;
-    static_assert(GlobalPosition::dimension == 3, "Only implemented for world dimension 3");
 
-    std::vector<GlobalPosition> points(numPoints);
+    // resize the points vector
+    const std::size_t numPoints = sincos.size()/2;
+    points.resize(numPoints);
 
     // make sure n is a unit vector
     auto n = normal;
     n /= n.two_norm();
 
-    // caculate a vector u perpendicular to n
+    // calculate a vector u perpendicular to n
     GlobalPosition u;
     if (abs(n[0]) < eps && abs(n[1]) < eps)
         if (abs(n[2]) < eps)
@@ -75,26 +83,47 @@ std::vector<GlobalPosition> circlePoints(const GlobalPosition& center,
     auto tangent = crossProduct(u, n);
     tangent *= radius/tangent.two_norm();
 
-    // the parameter with an offset
-    ctype t = 0 + 0.1;
     // insert the vertices
     for (std::size_t i = 0; i < numPoints; ++i)
     {
-        points[i] = GlobalPosition({u[0]*cos(t) + tangent[0]*sin(t) + center[0],
-                                    u[1]*cos(t) + tangent[1]*sin(t) + center[1],
-                                    u[2]*cos(t) + tangent[2]*sin(t) + center[2]});
+        points[i] = GlobalPosition({u[0]*sincos[2*i+1] + tangent[0]*sincos[2*i] + center[0],
+                                    u[1]*sincos[2*i+1] + tangent[1]*sincos[2*i] + center[1],
+                                    u[2]*sincos[2*i+1] + tangent[2]*sincos[2*i] + center[2]});
+    }
+}
 
-        t += 2*M_PI/numPoints;
+/*!
+ * \ingroup EmbeddedCoupling
+ * \brief returns a vector of points on a circle
+ * \param center the circle center
+ * \param normal the normal to the circle plane
+ * \param radius the circle radius
+ * \param numPoints the number of points
+ */
+template<class GlobalPosition, class Scalar>
+std::vector<GlobalPosition> circlePoints(const GlobalPosition& center,
+                                         const GlobalPosition& normal,
+                                         const Scalar radius,
+                                         const std::size_t numPoints = 20)
+{
+    std::vector<GlobalPosition> points;
 
-        // periodic t
-        if(t > 2*M_PI)
-            t -= 2*M_PI;
+    // precompute the sin/cos
+    std::vector<Scalar> sincos(2*numPoints);
+    Scalar t = 0 + 0.1; // add arbitrary offset
+    for (std::size_t i = 0; i < numPoints; ++i)
+    {
+        using std::sin; using std::cos;
+        sincos[2*i] = sin(t);
+        sincos[2*i + 1] = cos(t);
+        t += 2*M_PI/numPoints;
+        if(t > 2*M_PI) t -= 2*M_PI;
     }
 
+    circlePoints(points, sincos, center, normal, radius);
     return points;
 }
 
-} // end namespace EmbeddedCoupling
-} // end namespace Dumux
+} // end namespace Dumux::EmbeddedCoupling
 
 #endif
-- 
GitLab


From 2ddc53badaf05afee065246b5b2fea85c0d5d297 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Fri, 12 Jun 2020 10:05:33 +0200
Subject: [PATCH 2/2] [test] Add test for circlePoints helper function

---
 test/common/geometry/CMakeLists.txt       |  1 +
 test/common/geometry/test_circlepoints.cc | 44 +++++++++++++++++++++++
 2 files changed, 45 insertions(+)
 create mode 100644 test/common/geometry/test_circlepoints.cc

diff --git a/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt
index 90a842d59a..92ea708139 100644
--- a/test/common/geometry/CMakeLists.txt
+++ b/test/common/geometry/CMakeLists.txt
@@ -8,6 +8,7 @@ dumux_add_test(SOURCES test_2d2d_intersection.cc LABELS unit)
 dumux_add_test(SOURCES test_2d3d_intersection.cc LABELS unit)
 dumux_add_test(SOURCES test_graham_convex_hull.cc LABELS unit)
 dumux_add_test(SOURCES test_intersectingentity_cartesiangrid.cc LABELS unit)
+dumux_add_test(SOURCES test_circlepoints.cc LABELS unit)
 
 dune_symlink_to_source_files(FILES ball.msh)
 dumux_add_test(SOURCES test_intersectionentityset.cc
diff --git a/test/common/geometry/test_circlepoints.cc b/test/common/geometry/test_circlepoints.cc
new file mode 100644
index 0000000000..3b5a02750c
--- /dev/null
+++ b/test/common/geometry/test_circlepoints.cc
@@ -0,0 +1,44 @@
+#include <config.h>
+
+#include <cmath>
+#include <vector>
+#include <algorithm>
+#include <dune/common/exceptions.hh>
+#include <dune/common/float_cmp.hh>
+#include <dune/common/fvector.hh>
+
+#include <dumux/multidomain/embedded/circlepoints.hh>
+
+int main (int argc, char *argv[])
+{
+    using namespace Dumux;
+
+    using Point = Dune::FieldVector<double, 3>;
+    const Point center(0.0);
+    const Point normal({1.0, 1.0, 0.0}); // doesn't need to be a unit normal vector
+    const double radius = 3.0;
+    static constexpr std::size_t numPoints = 100;
+
+    // make a circle and and ellipse with equal major radii -> should yield the same points
+    const auto circle = EmbeddedCoupling::circlePoints(center, normal, radius, numPoints);
+    const auto checkRadius = [&](const auto& p){ return Dune::FloatCmp::eq(p.two_norm(), radius, 1e-7); };
+    if (!std::all_of(circle.begin(), circle.end(), checkRadius))
+        DUNE_THROW(Dune::Exception, "Not all point have the same radius!");
+
+    Point unitNormal = normal;
+    unitNormal /= unitNormal.two_norm();
+    bool correctOrientation = true;
+    for (int i = 1; i < circle.size(); ++i)
+    {
+        auto n = crossProduct(circle[i]-center, circle[i-1]-center);
+        n /= n.two_norm();
+        correctOrientation &= Dune::FloatCmp::eq(std::abs(n*unitNormal), 1.0, 1e-10);
+    }
+    if (!correctOrientation)
+        DUNE_THROW(Dune::Exception, "Not all points have the correct orientation!");
+
+    if (circle.size() != numPoints)
+        DUNE_THROW(Dune::Exception, "Wrong number of points!");
+
+    return 0;
+}
-- 
GitLab