simple_partial reffactored

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

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;
         }

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:

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

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";
         }
     }

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;
     }