diff --git a/CMakeLists.txt b/CMakeLists.txt
index 170a698417dbd755dee6305615011516f5060eb2..ec32c249e5046947220cd51a011d427a2ec7aa3a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -12,23 +12,23 @@ find_package(Eigen3 REQUIRED)
# Library
add_subdirectory(source/numsimulation)
-add_subdirectory(source/potential)
add_subdirectory(source/potential_rectangular)
add_subdirectory(source/potential_warped)
+add_subdirectory(source/potential)
add_subdirectory(source/scalar)
add_subdirectory(source/scalar_v2)
+add_subdirectory(source/simple_partial)
#add_subdirectory(source/simple)
-#add_subdirectory(source/simple_partial)
# Demo
-add_subdirectory(demo/potential)
add_subdirectory(demo/potential_rectangular)
add_subdirectory(demo/potential_warped)
+add_subdirectory(demo/potential)
add_subdirectory(demo/scalar)
add_subdirectory(demo/scalar_v2)
+add_subdirectory(demo/simple_partial)
#add_subdirectory(demo/simple)
-#add_subdirectory(demo/simple_partial)
-add_custom_target(demo DEPENDS potential_rectangular_demo_all potential_warped_demo_all potential_demo_all scalar_demo_all scalar_v2_demo_all)
+add_custom_target(demo DEPENDS potential_rectangular_demo_all potential_warped_demo_all potential_demo_all scalar_demo_all scalar_v2_demo_all simple_partial_demo_all)
#add_custom_target(demo DEPENDS potential_rectangular_demo_all potential_warped_demo_all potential_demo_all scalar_demo_all scalar_v2_demo_all simple_partial_demo_all simple_demo_all)
# Automation
diff --git a/demo/simple_partial/main.cpp b/demo/simple_partial/main.cpp
index a6f95902ba2e29defaf7b5efef69839dae08849a..d33653854087c25d4bf90e15edfe04e339cae711 100644
--- a/demo/simple_partial/main.cpp
+++ b/demo/simple_partial/main.cpp
@@ -45,10 +45,10 @@ int _main(int argc, char **argv)
b = ((u[1] - u[0]) - c * (sqr(max[1]) - sqr(min[1]))) / (max[1] - min[1]);
a = u[0] - b * min[1] - c * sqr(min[1]);
boundaries.resize(4);
- boundaries[0] = new Line(Eigen::Vector2d(max[0], max[1]), Eigen::Vector2d(max[0], min[1]), false); //Right
- boundaries[1] = new Line(Eigen::Vector2d(max[0], min[1]), Eigen::Vector2d(min[0], min[1]), false); //Bottom
- boundaries[2] = new Line(Eigen::Vector2d(min[0], min[1]), Eigen::Vector2d(min[0], max[1]), false); //Left
- boundaries[3] = new Line(Eigen::Vector2d(min[0], max[1]), Eigen::Vector2d(max[0], max[1]), false); //Top
+ boundaries[0].figure = new Line(Eigen::Vector2d(max[0], max[1]), Eigen::Vector2d(max[0], min[1]), false); //Right
+ boundaries[1].figure = new Line(Eigen::Vector2d(max[0], min[1]), Eigen::Vector2d(min[0], min[1]), false); //Bottom
+ boundaries[2].figure = new Line(Eigen::Vector2d(min[0], min[1]), Eigen::Vector2d(min[0], max[1]), false); //Left
+ boundaries[3].figure = new Line(Eigen::Vector2d(min[0], max[1]), Eigen::Vector2d(max[0], max[1]), false); //Top
grid_origin = Eigen::Vector2d((min[0] + max[0]) / 2, (min[1] + max[1]) / 2);
exact = true;
}
@@ -68,7 +68,7 @@ int _main(int argc, char **argv)
w(reader.get_real("W"))
{
boundaries.resize(1);
- boundaries[0] = new Circle(Eigen::Vector2d::Zero(), reader.get_real("R"), true);
+ boundaries[0].figure = new Circle(Eigen::Vector2d::Zero(), reader.get_real("R"), true);
grid_origin = Eigen::Vector2d::Zero();
exact = true;
}
diff --git a/include/simple_partial/simple_partial.h b/include/simple_partial/simple_partial.h
index 8c83777c5048f35900a5e7d40f76e07a1a575a9e..43c2dfc36072609e3d2c0b4cf9b3412d2d5a4563 100644
--- a/include/simple_partial/simple_partial.h
+++ b/include/simple_partial/simple_partial.h
@@ -3,11 +3,12 @@ Developed as part of Numerical Flow Simulation course
Author: Kyrylo Sovailo */
#pragma once
-#include <vector>
-#include <unordered_set>
-#include <Eigen/Dense>
+
#include <numsimulation/common.h>
+#include <Eigen/Dense>
#include <fstream>
+#include <vector>
+#include <set>
using namespace numsimulation;
@@ -16,19 +17,22 @@ namespace simple_partial
struct Point;
struct Face;
struct Cell;
+ struct TemporaryCell;
/// Point. Cells and faces consist of points
struct Point
{
//Constants
Eigen::Vector2d coord; ///< Point coordinate
+
+ Point(Eigen::Vector2d coord); ///< Creates point
};
/// Face. Face is a line between two points
struct Face
{
/// Precalculated geometry data for each neighboring face
- struct Precalculated
+ struct FaceConstants
{
int neighbor_velocity_sign; ///< Sign of velocity in neighboring cell (1 or -1)
int neighbor_velocity_inwards_sign; ///< Sign of inwards velocity relative to velocity (1 or -1 or 0)
@@ -48,7 +52,7 @@ namespace simple_partial
double cell_cell_distance; ///< Precalculated distance between cells (corrected by normal)
double area; ///< Precalculated area
int pressure_sign; ///< Precalculated sign of pressure (1 or -1)
- std::vector<Precalculated> precalculated; ///< Precalculated geometry data
+ std::vector<FaceConstants> face_constants; ///< Precalculated geometry data
//Variables
double velocity; ///< Velocity
@@ -69,6 +73,8 @@ namespace simple_partial
//Variables
double pressure; ///< Pressure
+
+ Cell(double area, Eigen::Vector2d center); ///< Creates cell
};
/// Numerical method
@@ -81,11 +87,17 @@ namespace simple_partial
exponential
};
+ ///Boundary
+ struct Boundary
+ {
+ Figure *figure; ///< Boundary geometry
+ };
+
/// Parameters (parameter functions unimplemented)
struct AbstractParameters
{
// Geometry parameters
- std::vector<Boundary*> boundaries; ///< List of boundaries
+ std::vector<Boundary> boundaries; ///< List of boundaries
// Grid parameters
Eigen::Vector2d grid_origin;///< Grid origin, where the field probing begins
@@ -112,11 +124,24 @@ namespace simple_partial
const AbstractParameters *_parameters;
const std::string _filename;
std::ofstream _file;
- std::unordered_set<Point*> _points;
- std::unordered_set<Face*> _faces;
- std::unordered_set<Cell*> _cells;
+ std::set<Point*> _points;
+ std::set<Face*> _faces;
+ std::set<Cell*> _cells;
double _time = 0.0;
+
+ void _add_first_cell(std::map<Position, TemporaryCell> &active); ///< Adds first cell to active set
+ void _add_all_cells(std::map<Position, TemporaryCell> &active, std::map<Position, TemporaryCell> &passive); ///< Searches and adds cells to passive set
+ void _calculate_area(std::map<Position, TemporaryCell> &cells); ///< Calculate area of the cells
+ void _create_cells(std::map<Position, TemporaryCell> &cells); ///< Create cells
+ void _create_points(std::map<Position, TemporaryCell> &cells); ///< Create points
+ void _create_faces(std::map<Position, TemporaryCell> &cells); ///< Create faces
+ void _precalculate_faces(); ///< Precalculate face constants
+ void _interconnect_faces(std::map<Position, TemporaryCell> &cells); ///< Setup face->face pointers
+ void _precalculate_faces_cells(); ///< Precalculate face->cells constants
+ void _setup_velocity(); ///< Setup velocity on all faces
+ void _setup_fixed_velocity(); ///< Setup velocity on all fixed faces
+ void _setup_pressure(); ///< Setup pressure on all cells
double _a(double pe) const;
public:
diff --git a/source/scalar_v2/scalar_v2.cpp b/source/scalar_v2/scalar_v2.cpp
index d8d73f231b0000657e6c07efb1d1beede77ff063..f031f76144aca3cac4970508d17628e6b866eb39 100644
--- a/source/scalar_v2/scalar_v2.cpp
+++ b/source/scalar_v2/scalar_v2.cpp
@@ -6,7 +6,10 @@ Author: Kyrylo Sovailo */
#include <unordered_map>
/*
- Beginning with scalar_v2 all constructors are copypasted from scalar. Unique parts are marked with UNIQUE mark. Decent template/OOP refactoring is not the priority
+ Beginning with scalar_v2 all constructors are copypasted from scalar. Unique parts will be descried here and marked with UNIQUE mark
+
+ Unique:
+ 1. face_cell_distance is not calculated
*/
namespace scalar_v2
diff --git a/source/simple_partial/output.cpp b/source/simple_partial/output.cpp
index b8f43451c7b653ddd2b78b3632a6af98c4376a44..d8a67d53ffdc903167248c9eceee149d03b37f6a 100644
--- a/source/simple_partial/output.cpp
+++ b/source/simple_partial/output.cpp
@@ -9,29 +9,29 @@ void simple_partial::Solver::output()
_file << "T " << _time << "\n";
//Write pressure data
- for (std::unordered_set<Cell*>::iterator i = _cells.begin(); i != _cells.end(); i++)
+ for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
{
- Cell *p = *i;
- _file << "X " << p->center(0) << " Y " << p->center(1) << " EP " << _parameters->pressure(p->center, _time) << "\n";
+ Cell *cell = *icell;
+ _file << "X " << cell->center(0) << " Y " << cell->center(1) << " EP " << _parameters->pressure(cell->center, _time) << "\n";
}
//Write velocity data
- for (std::unordered_set<Face*>::iterator i = _faces.begin(); i != _faces.end(); i++)
+ for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
{
- Face *f = *i;
- _file << "X " << f->center(0) << " Y " << f->center(1);
- if (abs(f->normal(0)) > abs(f->normal(1)))
+ Face *face = *iface;
+ _file << "X " << face->center(0) << " Y " << face->center(1);
+ if (abs(face->normal(0)) > abs(face->normal(1)))
{
//Horizontal normal, vertical face
- _file << " VX " << (f->normal * f->velocity).x();
- if (_parameters->exact) _file << " EVX " << _parameters->velocity(f->center, _time).x();
+ _file << " VX " << (face->normal * face->velocity).x();
+ if (_parameters->exact) _file << " EVX " << _parameters->velocity(face->center, _time).x();
_file << "\n";
}
else
{
//Vertical normal, horizontal face
- _file << " VY " << (f->normal * f->velocity).y();
- if (_parameters->exact) _file << " EVY " << _parameters->velocity(f->center, _time).y();
+ _file << " VY " << (face->normal * face->velocity).y();
+ if (_parameters->exact) _file << " EVY " << _parameters->velocity(face->center, _time).y();
_file << "\n";
}
}
diff --git a/source/simple_partial/simple_partial.cpp b/source/simple_partial/simple_partial.cpp
index d8c6620c7076a21d35ad84eab32ecdd02a582a98..695310a9f9919d0fbf7638c571744d222844505a 100644
--- a/source/simple_partial/simple_partial.cpp
+++ b/source/simple_partial/simple_partial.cpp
@@ -5,6 +5,16 @@ Author: Kyrylo Sovailo */
#include <simple_partial/simple_partial.h>
#include <unordered_map>
+/*
+ Beginning with scalar_v2 all constructors are copypasted from simple_partial. Unique parts will be descried here and marked with UNIQUE mark
+
+ Unique:
+ 1. Faces have neighbors
+ 2. Faces have precalculated constants, cells do not
+ 3. There are no boundary conditions
+ 4. Only square grid
+*/
+
namespace simple_partial
{
enum class PointStatus
@@ -37,8 +47,9 @@ namespace simple_partial
{
std::vector<TemporaryPoint> points;
std::vector<TemporaryFace> faces;
- double area = 0;
- Eigen::Vector2d center = Eigen::Vector2d::Zero();
+ double area;
+ Eigen::Vector2d center;
+ bool fixed_scalar = false;
Cell *cell = nullptr;
TemporaryCell(GridType grid_type);
@@ -52,138 +63,138 @@ simple_partial::TemporaryCell::TemporaryCell(GridType grid_type)
else { points.resize(4); faces.resize(4); }
}
-simple_partial::Solver::Solver(const AbstractParameters *parameters, std::string filename) : _parameters(parameters), _filename(filename)
-{
- //Works like the scalar algorithm wit two differences:
- //1. Interconnects faces
- //2. Extended boundary conditions
+simple_partial::Point::Point(Eigen::Vector2d coord) : coord(coord) {}
- //Define sets, add first active point
- std::unordered_map<Position, TemporaryCell> passive, active;
+simple_partial::Cell::Cell(double area, Eigen::Vector2d center) : area(area), center(center) {}
+
+void simple_partial::Solver::_add_first_cell(std::map<Position, TemporaryCell> &active)
+{
+ TemporaryCell new_cell(GridType::square);
+ new_cell.points[0].status = PointStatus::active;
+ active.insert({Position(), new_cell});
+ const std::vector<PointNeighbor> neighbors = get_point_neighbors(Position(), GridType::square, 0);
+ for (unsigned int n = 0; n < neighbors.size(); n++)
{
- TemporaryCell new_cell(GridType::square);
- new_cell.points[0].status = PointStatus::active;
- active.insert({Position(), new_cell});
- std::vector<std::pair<Position, unsigned int>> neighbors = get_point_neighbors(Position(), GridType::square, 0);
- for (unsigned int n = 0; n < neighbors.size(); n++)
- {
- new_cell = TemporaryCell(GridType::square);
- new_cell.points[neighbors[n].second].status = PointStatus::active;
- active.insert({neighbors[n].first, new_cell});
- }
+ new_cell = TemporaryCell(GridType::square);
+ new_cell.points[neighbors[n].point].status = PointStatus::active;
+ active.insert({neighbors[n].position, new_cell});
}
+}
- //Search all other elements
+void simple_partial::Solver::_add_all_cells(std::map<Position, TemporaryCell> &active, std::map<Position, TemporaryCell> &passive)
+{
while (!active.empty())
{
//Iterate through active, probe unprobed faces, add new cells to to_be_active
- std::unordered_map<Position, TemporaryCell> to_be_active;
- for (std::unordered_map<Position, TemporaryCell>::iterator i = active.begin(); i != active.end(); i++)
+ std::map<Position, TemporaryCell> to_be_active;
+ for (std::map<Position, TemporaryCell>::iterator cell = active.begin(); cell != active.end(); cell++)
{
//For each point
- std::vector<Eigen::Vector2d> points = get_points(i->first, GridType::square, parameters->grid_origin, parameters->grid_size);
+ std::vector<Eigen::Vector2d> points = get_points(cell->first, GridType::square, _parameters->grid_origin, _parameters->grid_size);
for (unsigned int p = 0; p < get_shape(GridType::square); p++)
{
- if (i->second.points[p].status != PointStatus::active) continue;
+ if (cell->second.points[p].status != PointStatus::active) continue;
- //Trying to probe counterclockwise
- unsigned int next_ccw = ((p == (i->second.points.size() - 1)) ? 0 : (p + 1));
- if (!i->second.faces[p].probed)
+ //Trying to probe counterclockwise
+ const unsigned int next_ccw = ((p == (cell->second.points.size() - 1)) ? 0 : (p + 1));
+ if (!cell->second.faces[p].probed)
{
Intersection intersection;
- for (unsigned int b = 0; b < parameters->boundaries.size(); b++)
+ //UNIQUE: no fixed_scalar
+ for (std::vector<Boundary>::const_iterator boundary = _parameters->boundaries.begin(); boundary != _parameters->boundaries.end(); boundary++)
{
- Intersection new_intersection = parameters->boundaries[b]->intersection(points[p], points[next_ccw]);
+ Intersection new_intersection = boundary->figure->intersection(points[p], points[next_ccw]);
if (new_intersection.valid && (!intersection.valid || ((new_intersection.coord-points[p]).squaredNorm() < (intersection.coord-points[p]).squaredNorm())))
{
intersection = new_intersection;
}
}
- i->second.faces[p].probed = true;
- i->second.faces[p].valid = intersection.valid;
- i->second.faces[p].intersection = intersection.coord;
- std::pair<Position, unsigned int> neighbor = get_face_neighbor(i->first, GridType::square, p);
- std::unordered_map<Position, TemporaryCell>::iterator f = active.find(neighbor.first);
- f->second.faces[neighbor.second].probed = true;
- f->second.faces[neighbor.second].valid = intersection.valid;
- f->second.faces[neighbor.second].intersection = intersection.coord;
+ cell->second.faces[p].probed = true;
+ cell->second.faces[p].valid = intersection.valid;
+ cell->second.faces[p].intersection = intersection.coord;
+ const FaceNeighbor neighbor = get_face_neighbor(cell->first, GridType::square, p);
+ std::map<Position, TemporaryCell>::iterator find = active.find(neighbor.position);
+ find->second.faces[neighbor.face].probed = true;
+ find->second.faces[neighbor.face].valid = intersection.valid;
+ find->second.faces[neighbor.face].intersection = intersection.coord;
}
- if (!i->second.faces[p].valid && i->second.points[next_ccw].status == PointStatus::unreached)
+ if (!cell->second.faces[p].valid && cell->second.points[next_ccw].status == PointStatus::unreached)
{
- i->second.points[next_ccw].status = PointStatus::to_be_active;
- std::vector<std::pair<Position, unsigned int>> neighbors = get_point_neighbors(i->first, GridType::square, next_ccw);
- for (std::vector<std::pair<Position, unsigned int>>::iterator neighbor = neighbors.begin(); neighbor != neighbors.end(); neighbor++)
+ cell->second.points[next_ccw].status = PointStatus::to_be_active;
+ const std::vector<PointNeighbor> neighbors = get_point_neighbors(cell->first, GridType::square, next_ccw);
+ for (std::vector<PointNeighbor>::const_iterator neighbor = neighbors.begin(); neighbor != neighbors.end(); neighbor++)
{
- if (to_be_active.find(neighbor->first) != to_be_active.end()) //to_be_active points remain to_be_active
+ if (to_be_active.find(neighbor->position) != to_be_active.end()) //to_be_active points remain to_be_active
{
- to_be_active.find(neighbor->first)->second.points[neighbor->second].status = PointStatus::to_be_active;
+ to_be_active.find(neighbor->position)->second.points[neighbor->point].status = PointStatus::to_be_active;
}
- else if (active.find(neighbor->first) != active.end()) //Active points remain active
+ else if (active.find(neighbor->position) != active.end()) //Active points remain active
{
- active.find(neighbor->first)->second.points[neighbor->second].status = PointStatus::to_be_active;
+ active.find(neighbor->position)->second.points[neighbor->point].status = PointStatus::to_be_active;
}
- else if (passive.find(neighbor->first) != passive.end()) //Passive points become to_be_active
+ else if (passive.find(neighbor->position) != passive.end()) //Passive points become to_be_active
{
- passive.find(neighbor->first)->second.points[neighbor->second].status = PointStatus::to_be_active;
- to_be_active.insert(*passive.find(neighbor->first));
- passive.erase(neighbor->first);
+ passive.find(neighbor->position)->second.points[neighbor->point].status = PointStatus::to_be_active;
+ to_be_active.insert(*passive.find(neighbor->position));
+ passive.erase(neighbor->position);
}
else //Create to_be_active point
{
TemporaryCell new_cell(GridType::square);
- new_cell.points[neighbor->second].status = PointStatus::to_be_active;
- to_be_active.insert({neighbor->first, new_cell});
+ new_cell.points[neighbor->point].status = PointStatus::to_be_active;
+ to_be_active.insert({neighbor->position, new_cell});
}
}
}
//Trying to probe clockwise
- unsigned int next_cw = ((p == 0) ? (i->second.points.size() - 1) : (p - 1));
- if (!i->second.faces[next_cw].probed)
+ unsigned int next_cw = ((p == 0) ? (cell->second.points.size() - 1) : (p - 1));
+ if (!cell->second.faces[next_cw].probed)
{
Intersection intersection;
- for (unsigned int b = 0; b < parameters->boundaries.size(); b++)
+ //UNIQUE: no fixed_scalar
+ for (std::vector<Boundary>::const_iterator boundary = _parameters->boundaries.begin(); boundary != _parameters->boundaries.end(); boundary++)
{
- Intersection new_intersection = parameters->boundaries[b]->intersection(points[p], points[next_cw]);
+ Intersection new_intersection = boundary->figure->intersection(points[p], points[next_cw]);
if (new_intersection.valid && (!intersection.valid || ((new_intersection.coord-points[p]).squaredNorm() < (intersection.coord-points[p]).squaredNorm())))
{
intersection = new_intersection;
}
}
- i->second.faces[next_cw].probed = true;
- i->second.faces[next_cw].valid = intersection.valid;
- i->second.faces[next_cw].intersection = intersection.coord;
- std::pair<Position, unsigned int> neighbor = get_face_neighbor(i->first, GridType::square, next_cw);
- std::unordered_map<Position, TemporaryCell>::iterator f = active.find(neighbor.first);
- f->second.faces[neighbor.second].probed = true;
- f->second.faces[neighbor.second].valid = intersection.valid;
- f->second.faces[neighbor.second].intersection = intersection.coord;
+ cell->second.faces[next_cw].probed = true;
+ cell->second.faces[next_cw].valid = intersection.valid;
+ cell->second.faces[next_cw].intersection = intersection.coord;
+ const FaceNeighbor neighbor = get_face_neighbor(cell->first, GridType::square, next_cw);
+ std::map<Position, TemporaryCell>::iterator find = active.find(neighbor.position);
+ find->second.faces[neighbor.face].probed = true;
+ find->second.faces[neighbor.face].valid = intersection.valid;
+ find->second.faces[neighbor.face].intersection = intersection.coord;
}
- if (!i->second.faces[next_cw].valid && i->second.points[next_cw].status == PointStatus::unreached)
+ if (!cell->second.faces[next_cw].valid && cell->second.points[next_cw].status == PointStatus::unreached)
{
- i->second.points[next_cw].status = PointStatus::to_be_active;
- std::vector<std::pair<Position, unsigned int>> neighbors = get_point_neighbors(i->first, GridType::square, next_cw);
- for (std::vector<std::pair<Position, unsigned int>>::iterator neighbor = neighbors.begin(); neighbor != neighbors.end(); neighbor++)
+ cell->second.points[next_cw].status = PointStatus::to_be_active;
+ const std::vector<PointNeighbor> neighbors = get_point_neighbors(cell->first, GridType::square, next_cw);
+ for (std::vector<PointNeighbor>::const_iterator neighbor = neighbors.begin(); neighbor != neighbors.end(); neighbor++)
{
- if (to_be_active.find(neighbor->first) != to_be_active.end()) //to_be_active points remain to_be_active
+ if (to_be_active.find(neighbor->position) != to_be_active.end()) //to_be_active points remain to_be_active
{
- to_be_active.find(neighbor->first)->second.points[neighbor->second].status = PointStatus::to_be_active;
+ to_be_active.find(neighbor->position)->second.points[neighbor->point].status = PointStatus::to_be_active;
}
- else if (active.find(neighbor->first) != active.end()) //Active points remain active
+ else if (active.find(neighbor->position) != active.end()) //Active points remain active
{
- active.find(neighbor->first)->second.points[neighbor->second].status = PointStatus::to_be_active;
+ active.find(neighbor->position)->second.points[neighbor->point].status = PointStatus::to_be_active;
}
- else if (passive.find(neighbor->first) != passive.end()) //Passive points become to_be_active
+ else if (passive.find(neighbor->position) != passive.end()) //Passive points become to_be_active
{
- passive.find(neighbor->first)->second.points[neighbor->second].status = PointStatus::to_be_active;
- to_be_active.insert(*passive.find(neighbor->first));
- passive.erase(neighbor->first);
+ passive.find(neighbor->position)->second.points[neighbor->point].status = PointStatus::to_be_active;
+ to_be_active.insert(*passive.find(neighbor->position));
+ passive.erase(neighbor->position);
}
else //Create to_be_active point
{
TemporaryCell new_cell(GridType::square);
- new_cell.points[neighbor->second].status = PointStatus::to_be_active;
- to_be_active.insert({neighbor->first, new_cell});
+ new_cell.points[neighbor->point].status = PointStatus::to_be_active;
+ to_be_active.insert({neighbor->position, new_cell});
}
}
}
@@ -194,33 +205,35 @@ simple_partial::Solver::Solver(const AbstractParameters *parameters, std::string
active.insert(to_be_active.begin(), to_be_active.end());
//Make to_be_active points active, make active points passive, divide active into new_active and passive
- std::unordered_map<Position, TemporaryCell> new_active;
- for (std::unordered_map<Position, TemporaryCell>::iterator i = active.begin(); i != active.end(); i++)
+ std::map<Position, TemporaryCell> new_active;
+ for (std::map<Position, TemporaryCell>::iterator cell = active.begin(); cell != active.end(); cell++)
{
bool active = false;
- for (unsigned int p = 0; p < i->second.points.size(); p++)
+ for (unsigned int p = 0; p < cell->second.points.size(); p++)
{
- if (i->second.points[p].status == PointStatus::to_be_active)
+ if (cell->second.points[p].status == PointStatus::to_be_active)
{
- i->second.points[p].status = PointStatus::active;
+ cell->second.points[p].status = PointStatus::active;
active = true;
}
- else if (i->second.points[p].status == PointStatus::active)
+ else if (cell->second.points[p].status == PointStatus::active)
{
- i->second.points[p].status = PointStatus::passive;
+ cell->second.points[p].status = PointStatus::passive;
}
}
- if (active) new_active.insert(*i);
- else passive.insert(*i);
+ if (active) new_active.insert(*cell);
+ else passive.insert(*cell);
}
active = new_active;
}
+}
- //Calculating cells' area
- for (std::unordered_map<Position, TemporaryCell>::iterator i = passive.begin(); i != passive.end(); i++)
+void simple_partial::Solver::_calculate_area(std::map<Position, TemporaryCell> &cells)
+{
+ for (std::map<Position, TemporaryCell>::iterator i = cells.begin(); i != cells.end(); i++)
{
//Create list of points
- std::vector<Eigen::Vector2d> points = get_points(i->first, GridType::square, parameters->grid_origin, parameters->grid_size);
+ std::vector<Eigen::Vector2d> points = get_points(i->first, GridType::square, _parameters->grid_origin, _parameters->grid_size);
std::vector<Eigen::Vector2d> point_list;
for (unsigned int p = 0; p < get_shape(GridType::square); p++)
{
@@ -230,6 +243,8 @@ simple_partial::Solver::Solver(const AbstractParameters *parameters, std::string
}
//Calculate area
+ i->second.area = 0;
+ i->second.center = Eigen::Vector2d::Zero();
for (unsigned int p = 1; p < (point_list.size() - 1); p++)
{
double a = (point_list[p] - point_list[0]).norm();
@@ -243,96 +258,106 @@ simple_partial::Solver::Solver(const AbstractParameters *parameters, std::string
}
i->second.center = i->second.center / i->second.area;
}
+}
+void simple_partial::Solver::_create_cells(std::map<Position, TemporaryCell> &cells)
+{
//Creating cells
- for (std::unordered_map<Position, TemporaryCell>::iterator i = passive.begin(); i != passive.end(); i++)
+ for (std::map<Position, TemporaryCell>::iterator cell = cells.begin(); cell != cells.end(); cell++)
{
- if (i->second.area >= parameters->area_threshold)
+ if (cell->second.area >= _parameters->area_threshold)
{
- _cells.insert(i->second.cell = new Cell);
- i->second.cell->pressure = parameters->pressure(i->second.center, 0);
- i->second.cell->center = i->second.center;
- i->second.cell->area = i->second.area;
+ _cells.insert(cell->second.cell = new Cell(cell->second.area, cell->second.center));
}
}
+}
+void simple_partial::Solver::_create_points(std::map<Position, TemporaryCell> &cells)
+{
//Creating points, faces and adding neighbors
- for (std::unordered_map<Position, TemporaryCell>::iterator i = passive.begin(); i != passive.end(); i++)
+ for (std::map<Position, TemporaryCell>::iterator cell = cells.begin(); cell != cells.end(); cell++)
{
- if (i->second.cell == nullptr) continue;
+ if (cell->second.cell == nullptr) continue;
//Creating points
- std::vector<Eigen::Vector2d> points = get_points(i->first, GridType::square, parameters->grid_origin, parameters->grid_size);
- for (unsigned int p = 0; p < i->second.points.size(); p++)
+ std::vector<Eigen::Vector2d> points = get_points(cell->first, GridType::square, _parameters->grid_origin, _parameters->grid_size);
+ for (unsigned int p = 0; p < cell->second.points.size(); p++)
{
//Regular points
- if (i->second.points[p].status == PointStatus::passive)
+ if (cell->second.points[p].status == PointStatus::passive)
{
- if (i->second.points[p].point == nullptr)
+ if (cell->second.points[p].point == nullptr)
{
- _points.insert(i->second.points[p].point = new Point());
- i->second.points[p].point->coord = points[p];
- std::vector<std::pair<Position, unsigned int>> neighbors = get_point_neighbors(i->first, GridType::square, p);
- for (std::vector<std::pair<Position, unsigned int>>::iterator neighbor = neighbors.begin(); neighbor != neighbors.end(); neighbor++)
- passive.find(neighbor->first)->second.points[neighbor->second].point = i->second.points[p].point;
+ _points.insert(cell->second.points[p].point = new Point(points[p]));
+ const std::vector<PointNeighbor> neighbors = get_point_neighbors(cell->first, GridType::square, p);
+ for (std::vector<PointNeighbor>::const_iterator neighbor = neighbors.begin(); neighbor != neighbors.end(); neighbor++)
+ cells.find(neighbor->position)->second.points[neighbor->point].point = cell->second.points[p].point;
}
- i->second.cell->points.push_back(i->second.points[p].point);
+ cell->second.cell->points.push_back(cell->second.points[p].point);
}
//Points on faces
- unsigned int next_ccw = ((p == (i->second.points.size() - 1)) ? 0 : (p + 1));
- if (i->second.points[p].status != i->second.points[next_ccw].status)
+ unsigned int next_ccw = ((p == (cell->second.points.size() - 1)) ? 0 : (p + 1));
+ if (cell->second.points[p].status != cell->second.points[next_ccw].status)
{
- if (i->second.faces[p].point == nullptr)
+ if (cell->second.faces[p].point == nullptr)
{
- _points.insert(i->second.faces[p].point = new Point());
- i->second.faces[p].point->coord = i->second.faces[p].intersection;
- std::pair<Position, unsigned int> neighbor = get_face_neighbor(i->first, GridType::square, p);
- passive.find(neighbor.first)->second.faces[neighbor.second].point = i->second.faces[p].point;
+ _points.insert(cell->second.faces[p].point = new Point(cell->second.faces[p].intersection));
+ const FaceNeighbor neighbor = get_face_neighbor(cell->first, GridType::square, p);
+ cells.find(neighbor.position)->second.faces[neighbor.face].point = cell->second.faces[p].point;
}
- i->second.cell->points.push_back(i->second.faces[p].point);
+ cell->second.cell->points.push_back(cell->second.faces[p].point);
}
}
+ }
+}
+
+void simple_partial::Solver::_create_faces(std::map<Position, TemporaryCell> &cells)
+{
+ //Creating points, faces and adding neighbors
+ for (std::map<Position, TemporaryCell>::iterator cell = cells.begin(); cell != cells.end(); cell++)
+ {
+ if (cell->second.cell == nullptr) continue;
//Creating faces
Face *irregular = nullptr;
- for (unsigned int p = 0; p < i->second.faces.size(); p++)
+ for (unsigned int p = 0; p < cell->second.faces.size(); p++)
{
- unsigned int next_ccw = ((p == (i->second.points.size() - 1)) ? 0 : (p + 1));
- if (i->second.points[p].status == PointStatus::passive || i->second.points[next_ccw].status == PointStatus::passive)
+ unsigned int next_ccw = ((p == (cell->second.points.size() - 1)) ? 0 : (p + 1));
+ if (cell->second.points[p].status == PointStatus::passive || cell->second.points[next_ccw].status == PointStatus::passive)
{
- std::pair<Position, unsigned int> neighbor = get_face_neighbor(i->first, GridType::square, p);
- if (i->second.faces[p].face == nullptr)
+ const FaceNeighbor neighbor = get_face_neighbor(cell->first, GridType::square, p);
+ if (cell->second.faces[p].face == nullptr)
{
- _faces.insert(i->second.faces[p].face = new Face());
- passive.find(neighbor.first)->second.faces[neighbor.second].face = i->second.faces[p].face;
+ _faces.insert(cell->second.faces[p].face = new Face());
+ cells.find(neighbor.position)->second.faces[neighbor.face].face = cell->second.faces[p].face;
}
- i->second.cell->faces.push_back(i->second.faces[p].face);
- i->second.cell->neighbors.push_back(passive.find(neighbor.first)->second.cell);
+ cell->second.cell->faces.push_back(cell->second.faces[p].face);
+ cell->second.cell->neighbors.push_back(cells.find(neighbor.position)->second.cell);
}
- if (i->second.points[p].status == PointStatus::passive && i->second.points[next_ccw].status == PointStatus::passive)
+ if (cell->second.points[p].status == PointStatus::passive && cell->second.points[next_ccw].status == PointStatus::passive)
{
- i->second.faces[p].face->point[0] = i->second.points[p].point;
- i->second.faces[p].face->point[1] = i->second.points[next_ccw].point;
+ cell->second.faces[p].face->point[0] = cell->second.points[p].point;
+ cell->second.faces[p].face->point[1] = cell->second.points[next_ccw].point;
}
- if (i->second.points[p].status == PointStatus::passive && i->second.points[next_ccw].status != PointStatus::passive)
+ if (cell->second.points[p].status == PointStatus::passive && cell->second.points[next_ccw].status != PointStatus::passive)
{
- i->second.faces[p].face->point[0] = i->second.points[p].point;
- i->second.faces[p].face->point[1] = i->second.faces[p].point;
+ cell->second.faces[p].face->point[0] = cell->second.points[p].point;
+ cell->second.faces[p].face->point[1] = cell->second.faces[p].point;
//Opening irregular face
_faces.insert(irregular = new Face);
- i->second.cell->faces.push_back(irregular);
- i->second.cell->neighbors.push_back(nullptr);
- irregular->point[0] = i->second.faces[p].point;
+ cell->second.cell->faces.push_back(irregular);
+ cell->second.cell->neighbors.push_back(nullptr);
+ irregular->point[0] = cell->second.faces[p].point;
}
- if (i->second.points[p].status != PointStatus::passive && i->second.points[next_ccw].status == PointStatus::passive)
+ if (cell->second.points[p].status != PointStatus::passive && cell->second.points[next_ccw].status == PointStatus::passive)
{
- i->second.faces[p].face->point[0] = i->second.faces[p].point;
- i->second.faces[p].face->point[1] = i->second.points[next_ccw].point;
+ cell->second.faces[p].face->point[0] = cell->second.faces[p].point;
+ cell->second.faces[p].face->point[1] = cell->second.points[next_ccw].point;
//Closing irregular face
if (irregular != nullptr)
{
- irregular->point[1] = i->second.faces[p].point;
+ irregular->point[1] = cell->second.faces[p].point;
irregular = nullptr;
}
}
@@ -340,64 +365,68 @@ simple_partial::Solver::Solver(const AbstractParameters *parameters, std::string
if (irregular != nullptr)
{
//Closing irregular face
- for (unsigned int p = 0; p < i->second.faces.size(); p++)
+ for (unsigned int p = 0; p < cell->second.faces.size(); p++)
{
- if (i->second.faces[p].point != nullptr) { irregular->point[1] = i->second.faces[p].point; break; }
+ if (cell->second.faces[p].point != nullptr) { irregular->point[1] = cell->second.faces[p].point; break; }
}
}
}
+}
- //Precalculating face-specific and setting velocity
- for (std::unordered_set<Face*>::iterator i = _faces.begin(); i != _faces.end(); i++)
+void simple_partial::Solver::_precalculate_faces()
+{
+ for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
{
- Face *face = *i;
- face->center = (face->point[0]->coord + face->point[1]->coord) / 2;
- Eigen::Vector2d face_vector = face->point[0]->coord - face->point[1]->coord; //should point same as normal
- face->length = face_vector.norm();
+ Face *face = *iface;
+ Eigen::Vector2d face_vector = face->point[0]->coord - face->point[1]->coord;
face->normal = rotate_cw(face_vector) / face_vector.norm();
-
- face->velocity = _parameters->velocity(face->center, 0).dot(face->normal);
+ face->length = face_vector.norm();
+ face->center = (face->point[0]->coord + face->point[1]->coord) / 2;
}
+}
- //Making reverse links
- for (std::unordered_map<Position, TemporaryCell>::iterator i = passive.begin(); i != passive.end(); i++)
+void simple_partial::Solver::_interconnect_faces(std::map<Position, TemporaryCell> &cells)
+{
+ for (std::map<Position, TemporaryCell>::iterator cell = cells.begin(); cell != cells.end(); cell++)
{
- if (i->second.cell == nullptr) continue;
+ if (cell->second.cell == nullptr) continue;
//Face to cell reverse pointer
- for (unsigned int f = 0; f < i->second.cell->faces.size(); f++) i->second.cell->faces[f]->cells.push_back(i->second.cell);
+ for (unsigned int f = 0; f < cell->second.cell->faces.size(); f++) cell->second.cell->faces[f]->cells.push_back(cell->second.cell);
//Face to face reverse pointer
- for (unsigned int f = 0; f < i->second.faces.size(); f++)
+ for (unsigned int f = 0; f < cell->second.faces.size(); f++)
{
- if (i->second.faces[f].face == nullptr) continue;
+ if (cell->second.faces[f].face == nullptr) continue;
- if (i->second.faces[(f+2)%4].face != nullptr) i->second.faces[f].face->neighbors.push_back(i->second.faces[(f+2)%4].face);
+ if (cell->second.faces[(f+2)%4].face != nullptr) cell->second.faces[f].face->neighbors.push_back(cell->second.faces[(f+2)%4].face);
- Position position[2] = { i->first, i->first };
+ Position position[2] = { cell->first, cell->first };
if (f == 0 || f == 2) { position[0].xi--; position[1].xi++; }
else { position[0].yi--; position[1].yi++; }
for (unsigned int j = 0; j < 2; j++)
{
- std::unordered_map<Position, TemporaryCell>::iterator found = passive.find(position[j]);
- if (found != passive.end() && found->second.faces[f].face != nullptr)
+ std::map<Position, TemporaryCell>::iterator find = cells.find(position[j]);
+ if (find != cells.end() && find->second.faces[f].face != nullptr)
{
bool match = false;
- for (unsigned int k = 0; k < i->second.faces[f].face->neighbors.size(); k++)
+ for (unsigned int k = 0; k < cell->second.faces[f].face->neighbors.size(); k++)
{
- if (i->second.faces[f].face->neighbors[k] == found->second.faces[f].face) match = true;
+ if (cell->second.faces[f].face->neighbors[k] == find->second.faces[f].face) match = true;
}
- if (!match) i->second.faces[f].face->neighbors.push_back(found->second.faces[f].face);
+ if (!match) cell->second.faces[f].face->neighbors.push_back(find->second.faces[f].face);
}
}
}
}
+}
- //Precalculating face- and neighbor-specific
- for (std::unordered_set<Face*>::iterator i = _faces.begin(); i != _faces.end(); i++)
+void simple_partial::Solver::_precalculate_faces_cells()
+{
+ for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
{
- Face *face = *i;
- if (face->cells.size() < 2) continue;
+ Face *face = *iface;
+ if (face->cells.size() < 2) continue; //Skip fixed
Eigen::Vector2d cell_cell_vector = face->cells[1]->center - face->cells[0]->center;
face->cell_cell_distance = cell_cell_vector.dot(face->normal);
@@ -405,59 +434,85 @@ simple_partial::Solver::Solver(const AbstractParameters *parameters, std::string
if (face->cell_cell_distance < 0) { cell_cell_vector *= -1; face->cell_cell_distance *= -1; face->pressure_sign *= -1; }
face->area = face->length * face->cell_cell_distance;
- face->precalculated.resize(face->neighbors.size());
+ face->face_constants.resize(face->neighbors.size());
for (unsigned int n = 0; n < face->neighbors.size(); n++)
{
- face->precalculated[n].neighbor_velocity_sign = (face->neighbors[n]->normal.dot(face->normal) < 0) ? (-1) : (1);
+ face->face_constants[n].neighbor_velocity_sign = (face->neighbors[n]->normal.dot(face->normal) < 0) ? (-1) : (1);
Eigen::Vector2d face_face_vector = face->center - face->neighbors[n]->center; //points inwards
double face_face_vector_y = face_face_vector.dot(face->normal);
double face_face_vector_x = face_face_vector.dot(rotate_cw(face->normal));
if (std::abs(face_face_vector_y) > std::abs(face_face_vector_x)) //upper or lower
{
- face->precalculated[n].neighbor_velocity_inwards_sign = ((face_face_vector_y < 0) ? (-1) : (1)) * face->precalculated[n].neighbor_velocity_sign; //negative if upper
- face->precalculated[n].face_face_distance = std::abs(face_face_vector_y);
- face->precalculated[n].border_length = face->length;
+ face->face_constants[n].neighbor_velocity_inwards_sign = ((face_face_vector_y < 0) ? (-1) : (1)) * face->face_constants[n].neighbor_velocity_sign; //negative if upper
+ face->face_constants[n].face_face_distance = std::abs(face_face_vector_y);
+ face->face_constants[n].border_length = face->length;
}
else //right or left
{
- face->precalculated[n].neighbor_velocity_inwards_sign = 0;
- face->precalculated[n].face_face_distance = std::abs(face_face_vector_x);
- face->precalculated[n].border_length = face->cell_cell_distance;
+ face->face_constants[n].neighbor_velocity_inwards_sign = 0;
+ face->face_constants[n].face_face_distance = std::abs(face_face_vector_x);
+ face->face_constants[n].border_length = face->cell_cell_distance;
}
}
}
+}
+
+void simple_partial::Solver::_setup_velocity()
+{
+ for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
+ {
+ Face *face = *iface;
+ face->velocity = _parameters->velocity(face->center, 0).dot(face->normal);
+ }
+}
+
+simple_partial::Solver::Solver(const AbstractParameters *parameters, std::string filename) : _parameters(parameters), _filename(filename)
+{
+ //Define sets, add first active point
+ std::map<Position, TemporaryCell> passive, active;
+ _add_first_cell(active);
+ _add_all_cells(active, passive);
+ _calculate_area(passive);
+ _create_cells(passive);
+ _create_points(passive);
+ _create_faces(passive);
+ _precalculate_faces();
+ _interconnect_faces(passive);
+ _precalculate_faces_cells();
+ _setup_velocity();
//Opening file
_file = std::ofstream(_filename, std::ios::binary);
- if (!_file.good()) throw std::runtime_error("scalar_v2::Solver::_output(): Cannot open file");
+ if (!_file.good()) throw std::runtime_error("simple_partial::Solver::_output(): Cannot open file");
//Write header
double min[2] = { std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity() };
double max[2] = { -std::numeric_limits<double>::infinity(), -std::numeric_limits<double>::infinity() };
- for (std::unordered_set<Cell*>::iterator i = _cells.begin(); i != _cells.end(); i++)
+ for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
{
- if ((*i)->center(0) < min[0]) min[0] = (*i)->center(0);
- if ((*i)->center(1) < min[1]) min[1] = (*i)->center(1);
- if ((*i)->center(0) > max[0]) max[0] = (*i)->center(0);
- if ((*i)->center(1) > max[1]) max[1] = (*i)->center(1);
+ Cell *cell = *icell;
+ if (cell->center(0) < min[0]) min[0] = cell->center(0);
+ if (cell->center(1) < min[1]) min[1] = cell->center(1);
+ if (cell->center(0) > max[0]) max[0] = cell->center(0);
+ if (cell->center(1) > max[1]) max[1] = cell->center(1);
}
_file << "NX " << (unsigned int)((max[0] - min[0]) / _parameters->grid_size[0]) << " NY " << (unsigned int)((max[1] - min[1]) / _parameters->grid_size[1]) << '\n';
}
simple_partial::Solver::~Solver()
{
- for (std::unordered_set<Point*>::iterator i = _points.begin(); i != _points.end(); i++)
+ for (std::set<Point*>::iterator i = _points.begin(); i != _points.end(); i++)
{
delete (*i);
}
- for (std::unordered_set<Face*>::iterator i = _faces.begin(); i != _faces.end(); i++)
+ for (std::set<Face*>::iterator i = _faces.begin(); i != _faces.end(); i++)
{
delete (*i);
}
- for (std::unordered_set<Cell*>::iterator i = _cells.begin(); i != _cells.end(); i++)
+ for (std::set<Cell*>::iterator i = _cells.begin(); i != _cells.end(); i++)
{
delete (*i);
}
diff --git a/source/simple_partial/solve.cpp b/source/simple_partial/solve.cpp
index a5484150e9e7db03914627eee80e9ec2fc472c0e..452376753b43c458f318f5cc913bbc905e7f6d2a 100644
--- a/source/simple_partial/solve.cpp
+++ b/source/simple_partial/solve.cpp
@@ -20,53 +20,51 @@ double simple_partial::Solver::_a(double pe) const
void simple_partial::Solver::step()
{
//Settting up pressure
- for (std::unordered_set<Cell*>::iterator i = _cells.begin(); i != _cells.end(); i++)
+ for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
{
- Cell *cell = *i;
+ Cell *cell = *icell;
cell->pressure = _parameters->pressure(cell->center, _time);
}
//Setting up fixed faces
- for (std::unordered_set<Face*>::iterator i = _faces.begin(); i != _faces.end(); i++)
+ for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
{
- Face *face = *i;
+ Face *face = *iface;
if (face->cells.size() < 2) face->velocity = _parameters->velocity(face->center, _time).dot(face->normal);
}
//Processing
- for (std::unordered_set<Face*>::iterator i = _faces.begin(); i != _faces.end(); i++)
+ for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
{
- Face *face = *i;
+ Face *face = *iface;
if (face->cells.size() < 2) continue; //Skip fixed cells
- double ap = _parameters->density * face->area / _parameters->step; //ap becomes ap_tilde till the end of the iteration
- double b = ap * face->velocity; //b becomes ap_tilde*next_velocity till the end of the iteration
+ double a0 = _parameters->density * face->area / _parameters->step; //a0 becomes a0_tilde till the end of the iteration
+ double b = a0 * face->velocity; //b becomes a0_tilde*next_velocity till the end of the iteration
b -= face->pressure_sign * (face->cells[1]->pressure - face->cells[0]->pressure) * face->area / face->cell_cell_distance;
for (unsigned int n = 0; n < face->neighbors.size(); n++)
{
//Precalculating
- //face velocity in the right direction, used as scalar
- double neighbor_velocity = face->precalculated[n].neighbor_velocity_sign * face->neighbors[n]->velocity;
- //face velocity in inwards direction, used as velocity
- double neighbor_velocity_inwards = face->precalculated[n].neighbor_velocity_inwards_sign * neighbor_velocity;
- double face_face_distance = face->precalculated[n].face_face_distance;
- double border_length = face->precalculated[n].border_length;
+ const double neighbor_velocity = face->face_constants[n].neighbor_velocity_sign * face->neighbors[n]->velocity;
+ const double neighbor_velocity_inwards = face->face_constants[n].neighbor_velocity_inwards_sign * neighbor_velocity;
+ const double face_face_distance = face->face_constants[n].face_face_distance;
+ const double border_length = face->face_constants[n].border_length;
//Kernel
- double d = _parameters->viscosity / face_face_distance;
- double pe = _parameters->density * neighbor_velocity_inwards * face_face_distance / _parameters->viscosity;
- double f = neighbor_velocity_inwards * _parameters->density;
- double a = d * border_length * _a(pe) + std::max(0.0, f * border_length);
- ap += a;
+ const double d = _parameters->viscosity / face_face_distance;
+ const double pe = _parameters->density * neighbor_velocity_inwards * face_face_distance / _parameters->viscosity;
+ const double f = neighbor_velocity_inwards * _parameters->density;
+ const double a = d * border_length * _a(pe) + std::max(0.0, f * border_length);
+ a0 += a;
b += a * neighbor_velocity;
}
- face->new_velocity = b / ap;
+ face->new_velocity = b / a0;
}
//Postprocessing
- for (std::unordered_set<Face*>::iterator i = _faces.begin(); i != _faces.end(); i++)
+ for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
{
- Face *face = *i;
+ Face *face = *iface;
if (face->cells.size() < 2) continue;
face->velocity = face->new_velocity;
}