Fixes

demo/simple/4.in

+SCENARIO PIPE
+
+GRID_SIZE_X 0.02 GRID_SIZE_Y 0.02
+AREA_THRESHOLD 0.000001
+LENGTH_THRESHOLD 0.000001
+METHOD LINEAR
+DENSITY 1000
+VISCOSITY 1
+STEP 0.01
+BETA 0.8
+MAX_ITERATION 1000
+MAX_E 0.0001
+P_MAX_RESIDUAL 0
+P_MAX_ITERATION 50
+P_SOLVER GAUSS
+P_REACTION IGNORE
+V_MAX_RESIDUAL 0
+V_MAX_ITERATION 50
+V_SOLVER GAUSS
+V_REACTION IGNORE
+
+STEPS 1000
+DIVIDER 500
\ No newline at end of file

demo/simple/CMakeLists.txt

@@ -26,5 +26,11 @@ add_custom_target(simple_demo3 DEPENDS ${PROJECT_SOURCE_DIR}/3.out)
 add_custom_command(OUTPUT ${PROJECT_SOURCE_DIR}/3.png COMMAND python3 ARGS ${CMAKE_SOURCE_DIR}/visualizer.py ${PROJECT_SOURCE_DIR}/3.out ${PROJECT_SOURCE_DIR}/3.png DEPENDS ${CMAKE_SOURCE_DIR}/visualizer.py ${PROJECT_SOURCE_DIR}/3.out VERBATIM)
 add_custom_target(simple_demo3v DEPENDS ${PROJECT_SOURCE_DIR}/3.png)
 
+# Demo 4
+add_custom_command(OUTPUT ${PROJECT_SOURCE_DIR}/4.out COMMAND $<TARGET_FILE:simple_demo> ${PROJECT_SOURCE_DIR}/4.in ${PROJECT_SOURCE_DIR}/4.out DEPENDS $<TARGET_FILE:simple_demo> ${PROJECT_SOURCE_DIR}/4.in VERBATIM)
+add_custom_target(simple_demo4 DEPENDS ${PROJECT_SOURCE_DIR}/4.out)
+add_custom_command(OUTPUT ${PROJECT_SOURCE_DIR}/4.png COMMAND python4 ARGS ${CMAKE_SOURCE_DIR}/visualizer.py ${PROJECT_SOURCE_DIR}/4.out ${PROJECT_SOURCE_DIR}/4.png DEPENDS ${CMAKE_SOURCE_DIR}/visualizer.py ${PROJECT_SOURCE_DIR}/4.out VERBATIM)
+add_custom_target(simple_demo4v DEPENDS ${PROJECT_SOURCE_DIR}/4.png)
+
 # Invokation
-add_custom_target(simple_demo_all DEPENDS simple_demo1v simple_demo2v simple_demo3v)
\ No newline at end of file
+add_custom_target(simple_demo_all DEPENDS simple_demo1v simple_demo2v simple_demo3v simple_demo4v)
\ No newline at end of file

demo/simple/main.cpp

@@ -73,13 +73,53 @@ int _main(int argc, char **argv)
         {
             boundaries.resize(1);
             boundaries[0].figure = new Circle(Eigen::Vector2d::Zero(), reader.get_real("R"), true);
+            boundaries[1].figure = new Line(Eigen::Vector2d(0,-0.5), Eigen::Vector2d(0,0.5), true);
             boundaries[0].fixed_pressure = false;
+            boundaries[1].fixed_pressure = true
             grid_origin = Eigen::Vector2d::Zero();
             exact = true;
         }
     };
+    struct PipeParameters : AbstractParameters
+    {
+        virtual Eigen::Vector2d velocity(Eigen::Vector2d coord, double time) const
+        {
+            return Eigen::Vector2d(0.1, 0);
+        }
+        virtual double pressure(Eigen::Vector2d coord, double time) const
+        {
+            return 0;
+        }
+        PipeParameters(ParameterReader &reader) : AbstractParameters(reader)
+        {
+            boundaries.resize(10);
+            boundaries[0].figure = new Line(Eigen::Vector2d(0, 1), Eigen::Vector2d(1.4, 1), false);
+            boundaries[1].figure = new Line(Eigen::Vector2d(1.4, 1), Eigen::Vector2d(1.4, 0.4), false);
+            boundaries[2].figure = new Line(Eigen::Vector2d(1.4, 0.4), Eigen::Vector2d(2, 0.4), false);
+            boundaries[3].figure = new Line(Eigen::Vector2d(2, 0.4), Eigen::Vector2d(2, 0), false);
+            boundaries[4].figure = new Line(Eigen::Vector2d(2, 0), Eigen::Vector2d(0.2, 0), false);
+            boundaries[5].figure = new Line(Eigen::Vector2d(0.2, 0), Eigen::Vector2d(0.2, 0.6), false);
+            boundaries[6].figure = new Line(Eigen::Vector2d(0.2, 0.6), Eigen::Vector2d(0.6, 0.6), false);
+            boundaries[7].figure = new Line(Eigen::Vector2d(0.6, 0.6), Eigen::Vector2d(0.6, 0.8), false);
+            boundaries[8].figure = new Line(Eigen::Vector2d(0.6, 0.8), Eigen::Vector2d(0, 0.8), false);
+            boundaries[9].figure = new Line(Eigen::Vector2d(0, 0.8), Eigen::Vector2d(0, 1), false);
+            boundaries[0].fixed_pressure = false;
+            boundaries[1].fixed_pressure = false;
+            boundaries[2].fixed_pressure = false;
+            boundaries[3].fixed_pressure = true;
+            boundaries[4].fixed_pressure = false;
+            boundaries[5].fixed_pressure = false;
+            boundaries[6].fixed_pressure = false;
+            boundaries[7].fixed_pressure = false;
+            boundaries[8].fixed_pressure = false;
+            boundaries[9].fixed_pressure = false;
+            grid_origin = Eigen::Vector2d(0.4, 0.4);
+            exact = false;
+        }
+    };
     if (reader.get_string("SCENARIO") == "POISEUILLE") parameters = new PoiseuilleParameters(reader);
     else if (reader.get_string("SCENARIO") == "ROTATION") parameters = new RotationParameters(reader);
+    else if (reader.get_string("SCENARIO") == "PIPE") parameters = new PipeParameters(reader);
     else throw std::runtime_error("_main(): Unrecognized scenario");
     
     //Read the rest of parameters

source/numsimulation/CMakeLists.txt

@@ -6,4 +6,4 @@ project(commmon)
 
 add_library(common SHARED common.cpp)
 target_link_libraries(common PUBLIC Eigen3::Eigen)
-target_include_directories(common PUBLIC ${CMAKE_SOURCE_DIR}/include)
\ No newline at end of file
+target_include_directories(common PUBLIC ${CMAKE_SOURCE_DIR}/imclude)
\ No newline at end of file

source/numsimulation/common.cpp

@@ -356,9 +356,9 @@ numsimulation::Line::Line(Eigen::Vector2d a, Eigen::Vector2d b, bool normal_cw)
 numsimulation::Intersection numsimulation::Line::intersection(Eigen::Vector2d a, Eigen::Vector2d b)
 {
     double t = ((a.x()-_a.x()) * (_a.y()-_b.y()) - (a.y()-_a.y()) * (_a.x()-_b.x())) / ((a.x()-b.x()) * (_a.y()-_b.y()) - (a.y()-b.y()) * (_a.x()-_b.x()));
-    if (t < 0 || t > 1) return Intersection();
+    if (!(t >= 0 && t <= 1)) return Intersection();
     double u = ((a.x()-_a.x()) * (a.y()-b.y()) - (a.y()-_a.y()) * (a.x()-b.x())) / ((a.x()-b.x()) * (_a.y()-_b.y()) - (a.y()-b.y()) * (_a.x()-_b.x()));
-    if (u < 0 || u > 1) return Intersection();
+    if (!(u >= 0 && u <= 1)) return Intersection();
     return Intersection(a * (1.0 - t) + b * t, _b - _a, _normal_cw ? rotate_cw(_b-_a) : rotate_ccw(_b-_a));
 }
 

source/simple/simple.cpp

@@ -317,7 +317,7 @@ void simple::Solver::_create_points(std::map<Position, TemporaryCell> &cells)
             {
                 if (cell->second.points[p].point == nullptr)
                 {
-                    _points.insert(cell->second.points[p].point = new Point(points[p]));
+                    _points.insert(cell->second.points[p].point);
                     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;