From ee4d2be67e0d791e10e3275bdc290f42c9b19b60 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Fri, 21 Nov 2014 15:06:37 +0100
Subject: [PATCH] Work in progress on AnisotropicEikonal.

---
 opm/core/tof/AnisotropicEikonal.cpp | 113 +++++++++++++++++++++++++++-
 opm/core/tof/AnisotropicEikonal.hpp |  11 ++-
 2 files changed, 121 insertions(+), 3 deletions(-)

diff --git a/opm/core/tof/AnisotropicEikonal.cpp b/opm/core/tof/AnisotropicEikonal.cpp
index d1eecc10f..784d8af3a 100644
--- a/opm/core/tof/AnisotropicEikonal.cpp
+++ b/opm/core/tof/AnisotropicEikonal.cpp
@@ -20,6 +20,7 @@
 #include <opm/core/tof/AnisotropicEikonal.hpp>
 #include <opm/core/grid/GridUtilities.hpp>
 #include <opm/core/grid.h>
+#include <set>
 
 namespace Opm
 {
@@ -34,6 +35,7 @@ namespace Opm
 	}
 	cell_neighbours_ = vertexNeighbours(grid);
 	orderCounterClockwise(grid, cell_neighbours_);
+	considered_.reserve(100);
     }
 
     /// Solve the eikonal equation.
@@ -42,12 +44,119 @@ namespace Opm
     /// \param[out] solution          Array of solution to the eikonal equation.
     void AnisotropicEikonal2d::solve(const double* metric,
 				     const std::vector<int>& startcells,
-				     std::vector<double>& solution) const
+				     std::vector<double>& solution)
     {
 	// The algorithm used is described in J.A. Sethian and A. Vladimirsky,
 	// "Ordered Upwind Methods for Static Hamilton-Jacobi Equations".
+	// Notation in comments is as used in that paper: U is the solution,
+	// and q is the boundary condition. One difference is that we talk about
+	// grid cells instead of mesh points.
+	//
+	// Algorithm summary:
+	// 1. Put all cells in Far. U_i = \inf.
+	// 2. Move the startcells to Accepted. U_i = q(x_i)
+	// 3. Move cells adjacent to startcells to Considered, evaluate
+	//    U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k}
+	// 4. Find the Considered cell with the smallest value: r.
+	// 5. Move cell r to Accepted. Update AcceptedFront.
+	// 6. Move cells adjacent to r from Far to Considered.
+	// 7. Recompute the value for all Considered cells within
+	//    distance h * F_2/F1 from x_r. Use min of previous and new.
+	// 8. If Considered is not empty, go to step 4.
+
+	// 1. Put all cells in Far. U_i = \inf.
+	const int num_cells = grid_.number_of_cells;
+	const double inf = 1e100;
+	solution.clear();
+	solution.resize(num_cells, inf);
+	considered_.clear();
+	is_considered_.clear();
+	is_considered_.resize(num_cells, false);
+
+	// 2. Move the startcells to Accepted. U_i = q(x_i)
+	std::vector<char> accepted(num_cells, false);
+	const int num_startcells = startcells.size();
+	for (int ii = 0; ii < num_startcells; ++ii) {
+	    accepted[startcells[ii]] = true;
+	    solution[startcells[ii]] = 0.0;
+	}
+	std::set<int> accepted_front(startcells.begin(), startcells.end());
+
+	// 3. Move cells adjacent to startcells to Considered, evaluate
+	//    U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k}
+	for (int ii = 0; ii < num_startcells; ++ii) {
+	    const int scell = startcells[ii];
+	    const int num_nb = cell_neighbours_[scell].size();
+	    for (int nb = 0; nb < num_nb; ++nb) {
+		const int nb_cell = cell_neighbours_[scell][nb];
+		if (!is_considered_[nb_cell]) {
+		    const double value = computeValue(nb_cell);
+		    pushConsidered(std::make_pair(value, nb_cell));
+		    is_considered_[nb_cell] = true;
+		}
+	    }
+	}
+
+	// 4. Find the Considered cell with the smallest value: r.
+
+	// 5. Move cell r to Accepted. Update AcceptedFront.
+
+	// 6. Move cells adjacent to r from Far to Considered.
+
+	// 7. Recompute the value for all Considered cells within
+	//    distance h * F_2/F1 from x_r. Use min of previous and new.
+
+	// 8. If Considered is not empty, go to step 4.
 
-	// 1. 
     }
 
+
+
+
+
+    double AnisotropicEikonal2d::computeValue(const int cell) const
+    {
+	const auto& nbs = cell_neighbours_[cell];
+	const int num_nbs = nbs.size();
+	double val = 1e100;
+	for (int ii = 0; ii < num_nbs; ++ii) {
+	    const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] };
+	    // if ... accepted front
+	}
+	return val;
+    }
+
+
+
+
+
+    const AnisotropicEikonal2d::ValueAndCell& AnisotropicEikonal2d::topConsidered() const
+    {
+	return considered_.front();
+    }
+
+
+
+
+
+    void AnisotropicEikonal2d::pushConsidered(const ValueAndCell& vc)
+    {
+	considered_.push_back(vc);
+	std::push_heap(considered_.begin(), considered_.end());
+    }
+
+
+
+
+
+    void AnisotropicEikonal2d::popConsidered()
+    {
+	std::pop_heap(considered_.begin(), considered_.end());
+	considered_.pop_back();
+    }
+
+
+
+
+
 } // namespace Opm
diff --git a/opm/core/tof/AnisotropicEikonal.hpp b/opm/core/tof/AnisotropicEikonal.hpp
index 1bf4b4025..6e6b73250 100644
--- a/opm/core/tof/AnisotropicEikonal.hpp
+++ b/opm/core/tof/AnisotropicEikonal.hpp
@@ -45,10 +45,19 @@ namespace Opm
         /// \param[out] solution          Array of solution to the eikonal equation.
         void solve(const double* metric,
 		   const std::vector<int>& startcells,
-		   std::vector<double>& solution) const;
+		   std::vector<double>& solution);
     private:
 	const UnstructuredGrid& grid_;
 	SparseTable<int> cell_neighbours_;
+	typedef std::pair<double, int> ValueAndCell;
+	std::vector<ValueAndCell> considered_;
+	std::vector<char> is_considered_;
+
+	double computeValue(const int cell) const;
+
+	const ValueAndCell& topConsidered() const;
+	void pushConsidered(const ValueAndCell& vc);
+	void popConsidered();
     };
 
 } // namespace Opm