Whitespace fix.

This commit is contained in:
Atgeirr Flø Rasmussen 2014-12-01 15:17:30 +01:00
parent 8cc68a231c
commit 2e2ad711e5
2 changed files with 88 additions and 88 deletions

View File

@ -28,13 +28,13 @@ namespace Opm
/// Construct solver. /// Construct solver.
/// \param[in] grid A 2d grid. /// \param[in] grid A 2d grid.
AnisotropicEikonal2d::AnisotropicEikonal2d(const UnstructuredGrid& grid) AnisotropicEikonal2d::AnisotropicEikonal2d(const UnstructuredGrid& grid)
: grid_(grid) : grid_(grid)
{ {
if (grid.dimensions != 2) { if (grid.dimensions != 2) {
OPM_THROW(std::logic_error, "Grid for AnisotropicEikonal2d must be 2d."); OPM_THROW(std::logic_error, "Grid for AnisotropicEikonal2d must be 2d.");
} }
cell_neighbours_ = cellNeighboursAcrossVertices(grid); cell_neighbours_ = cellNeighboursAcrossVertices(grid);
orderCounterClockwise(grid, cell_neighbours_); orderCounterClockwise(grid, cell_neighbours_);
} }
/// Solve the eikonal equation. /// Solve the eikonal equation.
@ -42,61 +42,61 @@ namespace Opm
/// \param[in] startcells Array of cells where u = 0 at the centroid. /// \param[in] startcells Array of cells where u = 0 at the centroid.
/// \param[out] solution Array of solution to the eikonal equation. /// \param[out] solution Array of solution to the eikonal equation.
void AnisotropicEikonal2d::solve(const double* metric, void AnisotropicEikonal2d::solve(const double* metric,
const std::vector<int>& startcells, const std::vector<int>& startcells,
std::vector<double>& solution) std::vector<double>& solution)
{ {
// The algorithm used is described in J.A. Sethian and A. Vladimirsky, // The algorithm used is described in J.A. Sethian and A. Vladimirsky,
// "Ordered Upwind Methods for Static Hamilton-Jacobi Equations". // "Ordered Upwind Methods for Static Hamilton-Jacobi Equations".
// Notation in comments is as used in that paper: U is the solution, // 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 // and q is the boundary condition. One difference is that we talk about
// grid cells instead of mesh points. // grid cells instead of mesh points.
// //
// Algorithm summary: // Algorithm summary:
// 1. Put all cells in Far. U_i = \inf. // 1. Put all cells in Far. U_i = \inf.
// 2. Move the startcells to Accepted. U_i = q(x_i) // 2. Move the startcells to Accepted. U_i = q(x_i)
// 3. Move cells adjacent to startcells to Considered, evaluate // 3. Move cells adjacent to startcells to Considered, evaluate
// U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k} // U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k}
// 4. Find the Considered cell with the smallest value: r. // 4. Find the Considered cell with the smallest value: r.
// 5. Move cell r to Accepted. Update AcceptedFront. // 5. Move cell r to Accepted. Update AcceptedFront.
// 6. Recompute the value for all Considered cells within // 6. Recompute the value for all Considered cells within
// distance h * F_2/F1 from x_r. Use min of previous and new. // distance h * F_2/F1 from x_r. Use min of previous and new.
// 7. Move cells adjacent to r from Far to Considered. // 7. Move cells adjacent to r from Far to Considered.
// 8. If Considered is not empty, go to step 4. // 8. If Considered is not empty, go to step 4.
// 1. Put all cells in Far. U_i = \inf. // 1. Put all cells in Far. U_i = \inf.
const int num_cells = grid_.number_of_cells; const int num_cells = grid_.number_of_cells;
const double inf = 1e100; const double inf = 1e100;
solution.clear(); solution.clear();
solution.resize(num_cells, inf); solution.resize(num_cells, inf);
is_accepted_.clear(); is_accepted_.clear();
is_accepted_.resize(num_cells, false); is_accepted_.resize(num_cells, false);
accepted_front_.clear(); accepted_front_.clear();
considered_.clear(); considered_.clear();
considered_handles_.clear(); considered_handles_.clear();
is_considered_.clear(); is_considered_.clear();
is_considered_.resize(num_cells, false); is_considered_.resize(num_cells, false);
// 2. Move the startcells to Accepted. U_i = q(x_i) // 2. Move the startcells to Accepted. U_i = q(x_i)
const int num_startcells = startcells.size(); const int num_startcells = startcells.size();
for (int ii = 0; ii < num_startcells; ++ii) { for (int ii = 0; ii < num_startcells; ++ii) {
is_accepted_[startcells[ii]] = true; is_accepted_[startcells[ii]] = true;
solution[startcells[ii]] = 0.0; solution[startcells[ii]] = 0.0;
} }
accepted_front_.insert(startcells.begin(), startcells.end()); accepted_front_.insert(startcells.begin(), startcells.end());
// 3. Move cells adjacent to startcells to Considered, evaluate // 3. Move cells adjacent to startcells to Considered, evaluate
// U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k} // U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k}
for (int ii = 0; ii < num_startcells; ++ii) { for (int ii = 0; ii < num_startcells; ++ii) {
const int scell = startcells[ii]; const int scell = startcells[ii];
const int num_nb = cell_neighbours_[scell].size(); const int num_nb = cell_neighbours_[scell].size();
for (int nb = 0; nb < num_nb; ++nb) { for (int nb = 0; nb < num_nb; ++nb) {
const int nb_cell = cell_neighbours_[scell][nb]; const int nb_cell = cell_neighbours_[scell][nb];
if (!is_accepted_[nb_cell] && !is_considered_[nb_cell]) { if (!is_accepted_[nb_cell] && !is_considered_[nb_cell]) {
const double value = computeValue(nb_cell, metric, solution.data()); const double value = computeValue(nb_cell, metric, solution.data());
pushConsidered(std::make_pair(value, nb_cell)); pushConsidered(std::make_pair(value, nb_cell));
} }
} }
} }
while (!considered_.empty()) { while (!considered_.empty()) {
// 4. Find the Considered cell with the smallest value: r. // 4. Find the Considered cell with the smallest value: r.
@ -178,17 +178,17 @@ namespace Opm
const double* solution) const const double* solution) const
{ {
// std::cout << "++++ computeValue(), cell = " << cell << std::endl; // std::cout << "++++ computeValue(), cell = " << cell << std::endl;
const auto& nbs = cell_neighbours_[cell]; const auto& nbs = cell_neighbours_[cell];
const int num_nbs = nbs.size(); const int num_nbs = nbs.size();
const double inf = 1e100; const double inf = 1e100;
double val = inf; double val = inf;
for (int ii = 0; ii < num_nbs; ++ii) { for (int ii = 0; ii < num_nbs; ++ii) {
const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] }; const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] };
if (accepted_front_.count(n[0]) && accepted_front_.count(n[1])) { if (accepted_front_.count(n[0]) && accepted_front_.count(n[1])) {
const double cand_val = computeFromTri(cell, n[0], n[1], metric, solution); const double cand_val = computeFromTri(cell, n[0], n[1], metric, solution);
val = std::min(val, cand_val); val = std::min(val, cand_val);
} }
} }
if (val == inf) { if (val == inf) {
// Failed to find two accepted front nodes adjacent to this, // Failed to find two accepted front nodes adjacent to this,
// so we go for a single-neighbour update. // so we go for a single-neighbour update.
@ -201,7 +201,7 @@ namespace Opm
} }
assert(val != inf); assert(val != inf);
// std::cout << "---> " << val << std::endl; // std::cout << "---> " << val << std::endl;
return val; return val;
} }
@ -214,18 +214,18 @@ namespace Opm
const int new_cell) const const int new_cell) const
{ {
// std::cout << "++++ computeValueUpdate(), cell = " << cell << std::endl; // std::cout << "++++ computeValueUpdate(), cell = " << cell << std::endl;
const auto& nbs = cell_neighbours_[cell]; const auto& nbs = cell_neighbours_[cell];
const int num_nbs = nbs.size(); const int num_nbs = nbs.size();
const double inf = 1e100; const double inf = 1e100;
double val = inf; double val = inf;
for (int ii = 0; ii < num_nbs; ++ii) { for (int ii = 0; ii < num_nbs; ++ii) {
const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] }; const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] };
if ((n[0] == new_cell || n[1] == new_cell) if ((n[0] == new_cell || n[1] == new_cell)
&& accepted_front_.count(n[0]) && accepted_front_.count(n[1])) { && accepted_front_.count(n[0]) && accepted_front_.count(n[1])) {
const double cand_val = computeFromTri(cell, n[0], n[1], metric, solution); const double cand_val = computeFromTri(cell, n[0], n[1], metric, solution);
val = std::min(val, cand_val); val = std::min(val, cand_val);
} }
} }
if (val == inf) { if (val == inf) {
// Failed to find two accepted front nodes adjacent to this, // Failed to find two accepted front nodes adjacent to this,
// so we go for a single-neighbour update. // so we go for a single-neighbour update.
@ -237,7 +237,7 @@ namespace Opm
} }
} }
// std::cout << "---> " << val << std::endl; // std::cout << "---> " << val << std::endl;
return val; return val;
} }
@ -336,7 +336,7 @@ namespace Opm
const AnisotropicEikonal2d::ValueAndCell& AnisotropicEikonal2d::topConsidered() const const AnisotropicEikonal2d::ValueAndCell& AnisotropicEikonal2d::topConsidered() const
{ {
return considered_.top(); return considered_.top();
} }
@ -345,7 +345,7 @@ namespace Opm
void AnisotropicEikonal2d::pushConsidered(const ValueAndCell& vc) void AnisotropicEikonal2d::pushConsidered(const ValueAndCell& vc)
{ {
HeapHandle h = considered_.push(vc); HeapHandle h = considered_.push(vc);
considered_handles_[vc.second] = h; considered_handles_[vc.second] = h;
is_considered_[vc.second] = true; is_considered_[vc.second] = true;
} }
@ -358,7 +358,7 @@ namespace Opm
{ {
is_considered_[considered_.top().second] = false; is_considered_[considered_.top().second] = false;
considered_handles_.erase(considered_.top().second); considered_handles_.erase(considered_.top().second);
considered_.pop(); considered_.pop();
} }

View File

@ -36,44 +36,44 @@ namespace Opm
class AnisotropicEikonal2d class AnisotropicEikonal2d
{ {
public: public:
/// Construct solver. /// Construct solver.
/// \param[in] grid A 2d grid. /// \param[in] grid A 2d grid.
explicit AnisotropicEikonal2d(const UnstructuredGrid& grid); explicit AnisotropicEikonal2d(const UnstructuredGrid& grid);
/// Solve the eikonal equation. /// Solve the eikonal equation.
/// \param[in] metric Array of metric tensors, M, for each cell. /// \param[in] metric Array of metric tensors, M, for each cell.
/// \param[in] startcells Array of cells where u = 0 at the centroid. /// \param[in] startcells Array of cells where u = 0 at the centroid.
/// \param[out] solution Array of solution to the eikonal equation. /// \param[out] solution Array of solution to the eikonal equation.
void solve(const double* metric, void solve(const double* metric,
const std::vector<int>& startcells, const std::vector<int>& startcells,
std::vector<double>& solution); std::vector<double>& solution);
private: private:
// Grid and topology. // Grid and topology.
const UnstructuredGrid& grid_; const UnstructuredGrid& grid_;
SparseTable<int> cell_neighbours_; SparseTable<int> cell_neighbours_;
// Keep track of accepted cells. // Keep track of accepted cells.
std::vector<char> is_accepted_; std::vector<char> is_accepted_;
std::set<int> accepted_front_; std::set<int> accepted_front_;
// Keep track of considered cells. // Keep track of considered cells.
typedef std::pair<double, int> ValueAndCell; typedef std::pair<double, int> ValueAndCell;
typedef boost::heap::compare<std::greater<ValueAndCell>> Comparator; typedef boost::heap::compare<std::greater<ValueAndCell>> Comparator;
typedef boost::heap::fibonacci_heap<ValueAndCell, Comparator> Heap; typedef boost::heap::fibonacci_heap<ValueAndCell, Comparator> Heap;
Heap considered_; Heap considered_;
typedef Heap::handle_type HeapHandle; typedef Heap::handle_type HeapHandle;
std::map<int, HeapHandle> considered_handles_; std::map<int, HeapHandle> considered_handles_;
std::vector<char> is_considered_; std::vector<char> is_considered_;
bool isClose(const int c1, const int c2, const double* metric) const; bool isClose(const int c1, const int c2, const double* metric) const;
double computeValue(const int cell, const double* metric, const double* solution) const; double computeValue(const int cell, const double* metric, const double* solution) const;
double computeValueUpdate(const int cell, const double* metric, const double* solution, const int new_cell) const; double computeValueUpdate(const int cell, const double* metric, const double* solution, const int new_cell) const;
double computeFromLine(const int cell, const int from, const double* metric, const double* solution) 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* solution) const; double computeFromTri(const int cell, const int n0, const int n1, const double* metric, const double* solution) const;
const ValueAndCell& topConsidered() const; const ValueAndCell& topConsidered() const;
void pushConsidered(const ValueAndCell& vc); void pushConsidered(const ValueAndCell& vc);
void popConsidered(); void popConsidered();
}; };
} // namespace Opm } // namespace Opm