From 8491d186d41ad40d4d609690a4a925ba301145e5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Wed, 3 Dec 2014 09:53:05 +0100
Subject: [PATCH] Make isClose() use actual grid and metric data.

Note that current implementation is somewhat ad-hoc, and not
in line with the algorithm of the paper.
---
 opm/core/tof/AnisotropicEikonal.cpp | 72 ++++++++++++++++++++++-------
 opm/core/tof/AnisotropicEikonal.hpp |  3 ++
 2 files changed, 58 insertions(+), 17 deletions(-)

diff --git a/opm/core/tof/AnisotropicEikonal.cpp b/opm/core/tof/AnisotropicEikonal.cpp
index f715e36a4..16c3f0feb 100644
--- a/opm/core/tof/AnisotropicEikonal.cpp
+++ b/opm/core/tof/AnisotropicEikonal.cpp
@@ -25,6 +25,37 @@
 namespace Opm
 {
 
+    namespace
+    {
+        /// Euclidean (isotropic) distance.
+        double distanceIso(const double v1[2],
+                           const double v2[2])
+        {
+            const double d[2] = { v2[0] - v1[0], v2[1] - v1[1] };
+            const double dist = std::sqrt(d[0]*d[0] + d[1]*d[1]);
+            return dist;
+        }
+
+        /// Anisotropic distance with respect to a metric g.
+        /// If d = v2 - v1, the distance is sqrt(d^T g d).
+        double distanceAniso(const double v1[2],
+                             const double v2[2],
+                             const double g[4])
+        {
+            const double d[2] = { v2[0] - v1[0], v2[1] - v1[1] };
+            const double dist = std::sqrt(+ g[0] * d[0] * d[0]
+                                          + g[1] * d[0] * d[1]
+                                          + g[2] * d[1] * d[0]
+                                          + g[3] * d[1] * d[1]);
+            return dist;
+        }
+    } // anonymous namespace
+
+
+
+
+
+
     /// Construct solver.
     /// \param[in] grid      A 2d grid.
     AnisotropicEikonal2d::AnisotropicEikonal2d(const UnstructuredGrid& grid)
@@ -35,6 +66,7 @@ namespace Opm
         }
         cell_neighbours_ = cellNeighboursAcrossVertices(grid);
         orderCounterClockwise(grid, cell_neighbours_);
+        computeGridRadius();
     }
 
     /// Solve the eikonal equation.
@@ -166,7 +198,10 @@ namespace Opm
                                        const int c2,
                                        const double* metric) const
     {
-        return true;
+        const double* v[] = { grid_.cell_centroids + 2*c1,
+                              grid_.cell_centroids + 2*c2 };
+        const double* m = metric + 4*c1;
+        return distanceAniso(v[0], v[1], m) < 3.0 * grid_radius_[c1];
     }
 
 
@@ -244,22 +279,6 @@ namespace Opm
 
 
 
-    double distanceAniso(const double v1[2],
-                         const double v2[2],
-                         const double g[4])
-    {
-        const double d[2] = { v2[0] - v1[0], v2[1] - v1[1] };
-        const double dist = std::sqrt(+ g[0] * d[0] * d[0]
-                                      + g[1] * d[0] * d[1]
-                                      + g[2] * d[1] * d[0]
-                                      + g[3] * d[1] * d[1]);
-        return dist;
-    }
-
-
-
-
-
     double AnisotropicEikonal2d::computeFromLine(const int cell,
                                                  const int from,
                                                  const double* metric,
@@ -364,5 +383,24 @@ namespace Opm
 
 
 
+    void AnisotropicEikonal2d::computeGridRadius()
+    {
+        const int num_cells = cell_neighbours_.size();
+        grid_radius_.resize(num_cells);
+        for (int cell = 0; cell < num_cells; ++cell) {
+            double radius = 0.0;
+            const double* v1 = grid_.cell_centroids + 2*cell;
+            const auto& nb = cell_neighbours_[cell];
+            for (auto it = nb.begin(); it != nb.end(); ++it) {
+                const double* v2 = grid_.cell_centroids + 2*(*it);
+                radius = std::max(radius, distanceIso(v1, v2));
+            }
+            grid_radius_[cell] = radius;
+        }
+    }
+
+
+
+
 
 } // namespace Opm
diff --git a/opm/core/tof/AnisotropicEikonal.hpp b/opm/core/tof/AnisotropicEikonal.hpp
index 215219e49..94daedb92 100644
--- a/opm/core/tof/AnisotropicEikonal.hpp
+++ b/opm/core/tof/AnisotropicEikonal.hpp
@@ -55,6 +55,7 @@ namespace Opm
         // Keep track of accepted cells.
         std::vector<char> is_accepted_;
         std::set<int> accepted_front_;
+        std::vector<double> grid_radius_;
 
         // Keep track of considered cells.
         typedef std::pair<double, int> ValueAndCell;
@@ -74,6 +75,8 @@ namespace Opm
         const ValueAndCell& topConsidered() const;
         void pushConsidered(const ValueAndCell& vc);
         void popConsidered();
+
+        void computeGridRadius();
     };
 
 } // namespace Opm