cleaning

classes/sph0/louis/src_cpp/Model.cpp

@@ -6,7 +6,7 @@
 #include <iomanip>
 #include "FixedParticle.h"
 #include "MobileParticle.h"
-#include "ParticleSorter.h"
+#include "Sorter.h"
 #include "Kernels.h"
 
 Model::Model() : sorter(*this)
@@ -74,22 +74,7 @@ Model::initialise()
     }
     file2.close();
 
-    // create kernel
-    switch (this->kernelKind)
-    {
-    case K_CUBIC_SPLINE:
-        this->kernel = new CubicSplineKernel();
-        break;
-    case K_QUADRATIC:
-        this->kernel = new QuadraticKernel();
-        break;
-    case K_QUINTIC_SPLINE:
-        this->kernel = new QuinticSplineKernel();
-        break;
-    default:
-        throw std::runtime_error("Bad value of kernel kind");
-    }
-    this->kappa = this->kernel->kappa;
+
 
     std::cout << "Initialisation finished." << std::endl;
     timers["initialisation"].stop();
@@ -190,7 +175,8 @@ Model::load_parameters(std::string const &param_path)
     file >> this->c_0;
     file >> this->rho_0;
     file >> this->dom_dim;
-    file >> this->kernelKind;
+    int kernelKind;
+    file >> kernelKind;
     file >> this->alpha;
     file >> this->beta;
     file >> this->eqnState;
@@ -200,6 +186,23 @@ Model::load_parameters(std::string const &param_path)
     file >> this->maxTime;
     file >> this->saveInt;
     file.close();
+
+    // create kernel
+    switch (kernelKind)
+    {
+    case K_CUBIC_SPLINE:
+        this->kernel = new CubicSplineKernel();
+        break;
+    case K_QUADRATIC:
+        this->kernel = new QuadraticKernel();
+        break;
+    case K_QUINTIC_SPLINE:
+        this->kernel = new QuinticSplineKernel();
+        break;
+    default:
+        throw std::runtime_error("Bad value of kernel kind");
+    }
+    this->kappa = this->kernel->kappa;    
 }
 
 /// Save a particle set onto disk.

classes/sph0/louis/src_cpp/Model.h

@@ -2,7 +2,7 @@
 #define SPH_MODEL_H
 
 #include "sph.h"
-#include "ParticleSorter.h"
+#include "Sorter.h"
 #include <vector>
 
 /// This class is used to manage all the particles,
@@ -13,15 +13,15 @@
 class Model
 {
 public:
-    ParticleSorter sorter; ///< sorter machine
+    Sorter sorter; ///< sorter machine
     Kernel *kernel;
 
     std::vector<Particle *> particles; ///< array of pointers toward particles
 
     int numFP;            ///< number of fixed particles
     int numMP;            ///< number of mobile particles
     int numPart;          ///< number of particles (FP+MP)
-    int kernelKind;       ///< kind of kernel
+
     int kappa;            ///< kappa linked to the eqn state
     double alpha;         ///< weighting factor in the artificial viscosity formulation
     double beta;          ///< weighting factor in the artificial viscosity formulation

classes/sph0/louis/src_cpp/Particle.cpp

@@ -102,13 +102,13 @@ void
 Particle::getNeighbours()
 {
     int RKstep = this->model.RKstep;
-    Eigen::Vector3d xyz = this->coord[RKstep]; // position of the particle
+    Eigen::Vector3d &xyz = this->coord[RKstep]; // position of the particle
 
     if (RKstep == 0)
     {
-        ParticleSorter *sorter = &this->model.sorter;
+        Sorter *sorter = &this->model.sorter;
         int nx = this->model.sorter.nx; // number of cells on a row of the domain
-        int cellsToCheck[27];                             // number of the cells to check for the neighbours
+        int cellsToCheck[27];           // number of the cells to check for the neighbours
 
         // calculates the number of the cell in which the particle is
         int xCell = (int)((xyz(0) - fmod(xyz(0), sorter->dx)) / sorter->dx) + 1;
@@ -169,14 +169,13 @@ Particle::getNeighbours()
                 std::vector<Particle *> *cells = &sorter->cells[cellsToCheck[i]];
                 for (size_t j = 0; j < cells->size(); j++)
                 {
-                    Particle *neigh = (*cells)[j];
-                    Eigen::Vector3d &neighXYZ = neigh->coord[RKstep];
+                    Particle *p = (*cells)[j];
+                    Eigen::Vector3d &neighXYZ = p->coord[RKstep];
                     double r = (xyz - neighXYZ).norm();
                     if (r <= this->model.kappa * this->h)
                     {
-                        // if (r > 1e-12) // Louis
-                        if (neigh != this)
-                            this->neighbours.push_back(Neighbour(neigh, r));
+                        if (p != this)
+                            this->neighbours.push_back(Neighbour(p, r));
                         else
                             twice++;
                     }
@@ -204,9 +203,9 @@ Particle::getNeighbours()
         // RK step2 - same neighbours and r is updated
         for (size_t i = 0; i < this->neighbours.size(); i++)
         {
-            Neighbour *link = &this->neighbours[i];
-            Eigen::Vector3d neighXYZ = link->p->coord[RKstep];
-            link->r = (xyz - neighXYZ).norm();
+            Neighbour *neigh = &this->neighbours[i];
+            Eigen::Vector3d &neighXYZ = neigh->p->coord[RKstep];
+            neigh->r = (xyz - neighXYZ).norm();
         }
     }
 }
@@ -230,97 +229,11 @@ Particle::gradW()
     {
         double r = this->neighbours[i].r;
         Particle *neigh = this->neighbours[i].p;
-        if(r>0.0)
-            this->vec_gradW[i] = (this->coord[RKstep] - neigh->coord[RKstep]) / r * this->model.kernel->dW(r, h);  
+        if (r > 0.0)
+            this->vec_gradW[i] = (this->coord[RKstep] - neigh->coord[RKstep]) / r * this->model.kernel->dW(r, h);
         else
             this->vec_gradW[i] = Eigen::Vector3d::Zero();
     }
-/*
-    switch (this->model.kernelKind)
-    {
-    case K_CUBIC_SPLINE:
-    {
-        double alpha_d = 3.0 / (2.0 * M_PI * h * h * h);
-        // values of alpha_d in table (2.1) p 23
-        for (size_t i = 0; i < this->neighbours.size(); i++)
-        {
-            double r = this->neighbours[i].r;
-            Particle *neigh = this->neighbours[i].p;
-            if ((r / h >= 0.0) && (r / h < 1.0))
-            {
-                this->vec_gradW[i] = alpha_d / h *
-                                     (1.5 * (r / h) * (r / h) - 2.0 * (r / h)) *
-                                     (this->coord[RKstep] - neigh->coord[RKstep]) / r;
-            }
-            else if ((r / h >= 1.0) && (r / h < 2.0))
-            {
-                this->vec_gradW[i] = alpha_d / h *
-                                     (-0.5 * (2.0 - r / h) * (2.0 - r / h)) *
-                                     (this->coord[RKstep] - neigh->coord[RKstep]) / r;
-            }
-            else
-            {
-                this->vec_gradW[i] = Eigen::Vector3d::Zero();
-            }
-        }
-        break;
-    }
-    case K_QUADRATIC:
-    {
-        double alpha_d = 5.0 / (4.0 * M_PI * h * h * h);
-        for (size_t i = 0; i < this->neighbours.size(); i++)
-        {
-            double r = this->neighbours[i].r;
-            Particle *neigh = this->neighbours[i].p;
-            if ((r / h >= 0.0) && (r / h <= 2.0))
-                this->vec_gradW[i] = alpha_d / h *
-                                     (0.375 * r / h - 0.75) *
-                                     (this->coord[RKstep] - neigh->coord[RKstep]) / r;
-
-            else
-                this->vec_gradW[i] = Eigen::Vector3d::Zero();
-        }
-        break;
-    }
-    case K_QUINTIC_SPLINE:
-    {
-        double alpha_d = 3.0 / (359.0 * M_PI * h * h * h);
-        for (size_t i = 0; i < this->neighbours.size(); i++)
-        {
-            double r = this->neighbours[i].r;
-            Particle *neigh = this->neighbours[i].p;
-            if ((r / h >= 0.0) && (r / h < 1.0))
-            {
-                this->vec_gradW[i] = alpha_d / h *
-                                     (-5.0 * pow(3.0 - r / h, 4) +
-                                      30.0 * pow(2.0 - r / h, 4) -
-                                      75.0 * pow(1.0 - r / h, 4)) *
-                                     (this->coord[RKstep] - neigh->coord[RKstep]) / r;
-            }
-            else if ((r / h >= 1.0) && (r / h < 2.0))
-            {
-                this->vec_gradW[i] = alpha_d / h *
-                                     (-5.0 * pow(3.0 - r / h, 4) +
-                                      30.0 * pow(2.0 - r / h, 4)) *
-                                     (this->coord[RKstep] - neigh->coord[RKstep]) / r;
-            }
-            else if ((r / h >= 2.0) && (r / h < 3.0))
-            {
-                this->vec_gradW[i] = alpha_d / h *
-                                     (-5.0 * pow(3.0 - r / h, 4)) *
-                                     (this->coord[RKstep] - neigh->coord[RKstep]) / r;
-            }
-            else
-            {
-                this->vec_gradW[i] = Eigen::Vector3d::Zero();
-            }
-        }
-        break;
-    }
-    default:
-        throw std::runtime_error("Bad value for kernel kind (1,2,3)");
-    }
-    */
 }
 
 /// Takes into account the fact that the kernel may be truncated.
@@ -332,7 +245,6 @@ Particle::kernel_corr()
     // if(this->vec_gradW_mod.size()<this->neighbours.size())
     //     this->vec_gradW_mod.resize(this->neighbours.size());
 
-
     int RKstep = this->model.RKstep;
 
     Eigen::Matrix3d M;

classes/sph0/louis/src_cpp/ParticleSorter.cpp → classes/sph0/louis/src_cpp/Sorter.cpp

@@ -1,19 +1,19 @@
-#include "ParticleSorter.h"
+#include "Sorter.h"
 #include "Model.h"
 #include "FixedParticle.h"
 #include "Neighbour.h"
 #include <fstream>
 #include <iostream>
 
-ParticleSorter::ParticleSorter(Model &m) : model(m)
+Sorter::Sorter(Model &m) : model(m)
 {
 }
 
 /// Put every particle in their corresponding cell.
 /// This will be useful in order to find the neighbours.
 
 void
-ParticleSorter::execute()
+Sorter::execute()
 {
     if (this->cells.size() == 0)
         this->init_cells();
@@ -58,7 +58,7 @@ ParticleSorter::execute()
 // reallocated at each iteration.
 
 void
-ParticleSorter::init_cells()
+Sorter::init_cells()
 {
     // calculates the necessary number of cells on a side
     double hmax = this->compute_hmax();
@@ -88,7 +88,7 @@ ParticleSorter::init_cells()
 /// This is useful when it is not constant over the particles.
 
 double 
-ParticleSorter::compute_hmax()
+Sorter::compute_hmax()
 {
     double hmax = 0.0;
     for(auto &p : this->model.particles)

classes/sph0/louis/src_cpp/ParticleSorter.h → classes/sph0/louis/src_cpp/Sorter.h

@@ -1,25 +1,26 @@
-#ifndef SPH_PARTICLESORTER_H
-#define SPH_PARTICLESORTER_H
+#ifndef SPH_SORTER_H
+#define SPH_SORTER_H
 
 #include "sph.h"
 #include "Neighbour.h"
 
 /// This class is able to sort the particles with the Linked-List method. 
 /// A grid is generated and the particles are sorted in each cell.
 
-class ParticleSorter
+class Sorter
 {
     Model &model;
     double h_max;    ///< maximum smoothing length
 
 public:
     double dx; ///< length of a side of a cube
-    int nx;  ///< number of cells on a row
+    int nx;    ///< number of cells on a row
+
+    /// vector of particles in each cell
+    std::vector<std::vector<Particle *>> cells; 
 
-    std::vector<std::vector<Particle *>> cells; ///< vector of lists that contain
-                                          /// the particles in a cell
 public:
-    ParticleSorter(Model &m);
+    Sorter(Model &m);
 
     void execute();
 
@@ -28,4 +29,4 @@ class ParticleSorter
     double compute_hmax();
 };
 
-#endif // SPH_PARTICLESORTER_H
+#endif // SPH_SORTER_H

classes/sph0/louis/src_cpp/sph.h

@@ -6,7 +6,7 @@ class Neighbour;
 class FixedParticle;
 class MobileParticle;
 class Model;
-class ParticleSorter;
+class Sorter;
 class Kernel;
 class CubicSplineKernel;
 class QuadraticKernel;