Initial versions of compute functions done.

Still has bugs.

opm/core/tof/AnisotropicEikonal.cpp

@@ -20,7 +20,7 @@
 #include <opm/core/tof/AnisotropicEikonal.hpp>
 #include <opm/core/grid/GridUtilities.hpp>
 #include <opm/core/grid.h>
-#include <set>
+#include <opm/core/utility/RootFinders.hpp>
 
 namespace Opm
 {
@@ -92,7 +92,7 @@ namespace Opm
 	    for (int nb = 0; nb < num_nb; ++nb) {
 		const int nb_cell = cell_neighbours_[scell][nb];
 		if (!is_accepted_[nb_cell] && !is_considered_[nb_cell]) {
-		    const double value = computeValue(nb_cell, metric);
+		    const double value = computeValue(nb_cell, metric, solution.data());
 		    pushConsidered(std::make_pair(value, nb_cell));
 		}
 	    }
@@ -130,7 +130,7 @@ namespace Opm
                 const int nb_cell = *it;
                 if (!is_accepted_[nb_cell] && !is_considered_[nb_cell]) {
                     assert(solution[nb_cell] == inf);
-                    const double value = computeValue(nb_cell, metric);
+                    const double value = computeValue(nb_cell, metric, solution.data());
                     pushConsidered(std::make_pair(value, nb_cell));
                 }
             }
@@ -140,7 +140,7 @@ namespace Opm
             for (auto it = considered_.begin(); it != considered_.end(); ++it) {
                 const int ccell = it->second;
                 if (isClose(rcell, ccell, metric)) {
-                    const double value = computeValue(ccell, metric);
+                    const double value = computeValue(ccell, metric, solution.data());
                     if (value < it->first) {
                         // Update value for considered cell.
                         // Note that as solution values decrease, their
@@ -173,7 +173,8 @@ namespace Opm
 
 
     double AnisotropicEikonal2d::computeValue(const int cell,
-                                              const double* metric) const
+                                              const double* metric,
+                                              const double* solution) const
     {
 	const auto& nbs = cell_neighbours_[cell];
 	const int num_nbs = nbs.size();
@@ -182,7 +183,7 @@ namespace Opm
 	for (int ii = 0; ii < num_nbs; ++ii) {
 	    const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] };
             if (accepted_front_.count(n[0]) && accepted_front_.count(n[1])) {
-                const double cand_val = computeFromTri(cell, n[0], n[1], metric);
+                const double cand_val = computeFromTri(cell, n[0], n[1], metric, solution);
                 val = std::min(val, cand_val);
             }
 	}
@@ -191,7 +192,7 @@ namespace Opm
             // so we go for a single-neighbour update.
             for (int ii = 0; ii < num_nbs; ++ii) {
                 if (accepted_front_.count(nbs[ii])) {
-                    const double cand_val = computeFromLine(cell, nbs[ii], metric);
+                    const double cand_val = computeFromLine(cell, nbs[ii], metric, solution);
                     val = std::min(val, cand_val);
                 }
             }
@@ -204,28 +205,85 @@ 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) const
+                                                 const double* metric,
+                                                 const double* solution) const
     {
         assert(!is_accepted_[cell]);
         assert(is_accepted_[from]);
-        return 2.0;
+        // Applying the first fundamental form to compute geodesic distance.
+        // Using the metric of 'cell', not 'from'.
+        const double dist = distanceAniso(grid_.cell_centroids + 2 * cell,
+                                          grid_.cell_centroids + 2 * from,
+                                          metric + 4 * cell);
+        return solution[from] + dist;
     }
 
 
 
 
 
+    struct DistanceDerivative
+    {
+        const double* x1;
+        const double* x2;
+        const double* x;
+        double u1;
+        double u2;
+        const double* g;
+        double operator()(const double theta) const
+        {
+            double a[2] = { x[0] - (1-theta)*x1[0] + theta*x2[0], x[1] - (1-theta)*x1[1] + theta*x2[1] };
+            double b[2] = { x1[0] - x2[0], x1[1] - x2[1] };
+            return u2 - u1 + 2*(a[0]*b[0]*g[0] + a[0]*b[1]*g[1] + a[1]*b[0]*g[2] + a[1]*b[1]*g[3]);
+        }
+    };
+
+
+
+
+
     double AnisotropicEikonal2d::computeFromTri(const int cell,
                                                 const int n0,
                                                 const int n1,
-                                                const double* metric) const
+                                                const double* metric,
+                                                const double* solution) const
     {
         assert(!is_accepted_[cell]);
         assert(is_accepted_[n0]);
         assert(is_accepted_[n1]);
-        return 1.0;
+        DistanceDerivative dd;
+        dd.x1 = grid_.cell_centroids + 2 * n0;
+        dd.x2 = grid_.cell_centroids + 2 * n1;
+        dd.x = grid_.cell_centroids + 2 * cell;
+        dd.u1 = solution[n0];
+        dd.u2 = solution[n1];
+        dd.g = metric + 4 * cell;
+        int iter = 0;
+        const double theta = RegulaFalsi<ContinueOnError>::solve(dd, 0.0, 1.0, 15, 1e-8, iter);
+        const double xt[2] = { (1-theta)*dd.x1[0] + theta*dd.x2[0],
+                               (1-theta)*dd.x1[1] + theta*dd.x2[1] };
+        const double d1 = distanceAniso(dd.x1, dd.x, dd.g) + solution[n0];
+        const double d2 = distanceAniso(dd.x2, dd.x, dd.g) + solution[n1];
+        const double dt = distanceAniso(xt, dd.x, dd.g) + (1-theta)*solution[n0] + theta*solution[n1];
+        return std::min(d1, std::min(d2, dt));
     }
 
 

opm/core/tof/AnisotropicEikonal.hpp

@@ -66,9 +66,9 @@ namespace Opm
 	std::vector<char> is_considered_;
 
         bool isClose(const int c1, const int c2, const double* metric) const;
-	double computeValue(const int cell, const double* metric) const;
+	double computeValue(const int cell, const double* metric, const double* solution) const;
-	double computeFromLine(const int cell, const int from, const double* metric) const;
+	double computeFromLine(const int cell, const int from, const double* metric, const double* solution) const;
-	double computeFromTri(const int cell, const int n0, const int n1, const double* metric) const;
+	double computeFromTri(const int cell, const int n0, const int n1, const double* metric, const double* solution) const;
 
 	const ValueAndCell& topConsidered() const;
 	void pushConsidered(const ValueAndCell& vc);