From 1d98fc2f20cb38f1bc0d0899201095e14d020cdb Mon Sep 17 00:00:00 2001 From: Romain Boman Date: Thu, 7 Mar 2024 17:31:36 +0100 Subject: [PATCH] cleaning --- classes/sph0/louis/src_cpp/Model.cpp | 39 +++--- classes/sph0/louis/src_cpp/Model.h | 6 +- classes/sph0/louis/src_cpp/Particle.cpp | 112 ++---------------- .../{ParticleSorter.cpp => Sorter.cpp} | 10 +- .../src_cpp/{ParticleSorter.h => Sorter.h} | 17 +-- classes/sph0/louis/src_cpp/sph.h | 2 +- 6 files changed, 51 insertions(+), 135 deletions(-) rename classes/sph0/louis/src_cpp/{ParticleSorter.cpp => Sorter.cpp} (94%) rename classes/sph0/louis/src_cpp/{ParticleSorter.h => Sorter.h} (54%) diff --git a/classes/sph0/louis/src_cpp/Model.cpp b/classes/sph0/louis/src_cpp/Model.cpp index 593d7ea3..e5fbcb5b 100644 --- a/classes/sph0/louis/src_cpp/Model.cpp +++ b/classes/sph0/louis/src_cpp/Model.cpp @@ -6,7 +6,7 @@ #include #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 ¶m_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 ¶m_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. diff --git a/classes/sph0/louis/src_cpp/Model.h b/classes/sph0/louis/src_cpp/Model.h index cdb2e065..6dfe2139 100644 --- a/classes/sph0/louis/src_cpp/Model.h +++ b/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 /// This class is used to manage all the particles, @@ -13,7 +13,7 @@ class Model { public: - ParticleSorter sorter; ///< sorter machine + Sorter sorter; ///< sorter machine Kernel *kernel; std::vector particles; ///< array of pointers toward particles @@ -21,7 +21,7 @@ class Model 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 diff --git a/classes/sph0/louis/src_cpp/Particle.cpp b/classes/sph0/louis/src_cpp/Particle.cpp index d7716b16..cc8edd9d 100644 --- a/classes/sph0/louis/src_cpp/Particle.cpp +++ b/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 *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()neighbours.size()) // this->vec_gradW_mod.resize(this->neighbours.size()); - int RKstep = this->model.RKstep; Eigen::Matrix3d M; diff --git a/classes/sph0/louis/src_cpp/ParticleSorter.cpp b/classes/sph0/louis/src_cpp/Sorter.cpp similarity index 94% rename from classes/sph0/louis/src_cpp/ParticleSorter.cpp rename to classes/sph0/louis/src_cpp/Sorter.cpp index a42eec52..1ed046de 100644 --- a/classes/sph0/louis/src_cpp/ParticleSorter.cpp +++ b/classes/sph0/louis/src_cpp/Sorter.cpp @@ -1,11 +1,11 @@ -#include "ParticleSorter.h" +#include "Sorter.h" #include "Model.h" #include "FixedParticle.h" #include "Neighbour.h" #include #include -ParticleSorter::ParticleSorter(Model &m) : model(m) +Sorter::Sorter(Model &m) : model(m) { } @@ -13,7 +13,7 @@ ParticleSorter::ParticleSorter(Model &m) : model(m) /// 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) diff --git a/classes/sph0/louis/src_cpp/ParticleSorter.h b/classes/sph0/louis/src_cpp/Sorter.h similarity index 54% rename from classes/sph0/louis/src_cpp/ParticleSorter.h rename to classes/sph0/louis/src_cpp/Sorter.h index 99458688..6bb01046 100644 --- a/classes/sph0/louis/src_cpp/ParticleSorter.h +++ b/classes/sph0/louis/src_cpp/Sorter.h @@ -1,5 +1,5 @@ -#ifndef SPH_PARTICLESORTER_H -#define SPH_PARTICLESORTER_H +#ifndef SPH_SORTER_H +#define SPH_SORTER_H #include "sph.h" #include "Neighbour.h" @@ -7,19 +7,20 @@ /// 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> cells; - std::vector> 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 diff --git a/classes/sph0/louis/src_cpp/sph.h b/classes/sph0/louis/src_cpp/sph.h index 46c8b551..e568209e 100644 --- a/classes/sph0/louis/src_cpp/sph.h +++ b/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;