Particle.hpp

Go to the documentation of this file.

//
// Author(s)    : - Jean-luc CHARLES I2M-DuMAS/ENSAM Talence France
//                  <jean-luc.charles@ensam.eu>
//                - Damien ANDRE     SPCTS/ENS Ceramique industrielle, Limoges France
//                  <damien.andre@unilim.fr>
//                - Jeremie GIRARDOT I2M-DuMAS/ENSAM Talence France
//                  <jeremie.girardot@ensam.eu>
//                - Cedric Hubert    LAMIH/UVHC Valenciennes France
//                  <cedric.hubert@univ-valenciennes.fr>
//                - Ivan IORDANOFF   I2M-DuMAS-MPI/ENSAM Talence France
//                  <ivan.iordanoff@ensam.eu>
//
// Copyright (C) 2008-2016 JL. CHARLES, D. ANDRE, I. IORDANOFF, J. GIRARDOT
//
//
//
//
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 
 
 
#ifndef GranOO_SPH_Particle_hpp
#define GranOO_SPH_Particle_hpp
 
 
 
#include "GranOO3/Physic/Node.hpp"
#include "GranOO3/SPH/Interaction.hpp"
 
 
namespace GranOO3
{
  namespace SPH
  {
 
    class Particle : public Physic::Node, public Core::PropClass<Particle>, public Core::Register<Particle> {
 
      GRANOO_CLASS(SPH::Particle, Physic::Node);
      GRANOO_CLASS_DEFAULT_COLOR(.8, .0, .1, 1.);
 
    public:
      Particle(const Geom::Point&, double radius);
      Particle();
      virtual ~Particle();
 
 
      GRANOO_ACCESS(h, double, _h);
      GRANOO_ACCESS(radius   , double, _radius   );
      GRANOO_ACCESS(density  , double, _density  );
      GRANOO_ACCESS(pressure , double, _pressure );
      GRANOO_ACCESS(viscosity, double, _viscosity);
 
      void   set_mass(double);
      double get_mass() const ;
 
      double get_bounding_radius() const;
      
      size_t neighbour_number() const;
      std::vector<Interaction>& get_interaction();
      void clear_interaction();
      void add_interaction(Particle& part);
 
      // UPDATE LOCAL PARAM
      template<class Kernel> void update_density();
      void update_pressure(const double gazStiffness, const double gazdensity);
      template<class Kernel> void apply_pressure_force();
      template<class Kernel> void apply_viscous_force();
 
      // USEFULL
      virtual std::string info() const;
      virtual void draw();
 
      // IO
      virtual std::ostream& write_ascii (std::ostream& out) const;
      virtual std::istream& read_ascii (std::istream& in);
 
    private:
      //BOOST SERIALIZATION
      friend class boost::serialization::access;
      template<class Archive> void serialize(Archive&, const unsigned int);
 
    private:
      double _mass;
      double _h;
      double _radius;
      double _density;
      double _pressure;
      double _viscosity;
      std::vector<Interaction> interaction_;
    };
 
    inline size_t
    Particle::neighbour_number() const {
      return interaction_.size();
    }
 
    inline std::vector<Interaction>&
    Particle::get_interaction() {
      return interaction_;
    }
 
    inline void
    Particle::clear_interaction() {
      interaction_.clear();
    }
 
    inline void
    Particle::add_interaction(Particle& part) {
      interaction_.emplace_back(&part, part.get_position() - get_position());
    }
 
    template<class Kernel> void
    Particle::update_density() {
      const Geom::Vector vecnull(0., 0., 0.);
      _density = 0.;
      _density += _mass * Kernel::k_def(_h, vecnull);
 
      for (auto& inter: interaction_) {
        const SPH::Particle* part2 = inter.neighbour;
        const Geom::Vector& vec =  inter.distance;
        const double hh = part2->get_h();
        _density += part2->get_mass()* Kernel::k_def(hh, vec);
      }
    }
 
    inline void
    Particle::update_pressure(const double gazStiffness, const double gazdensity) {
      _pressure = gazStiffness * (_density - gazdensity);
    }
 
    template<class Kernel> void
    Particle::apply_pressure_force() {
      const double pressure1 = _pressure;
      for (auto& inter: interaction_) {
        SPH::Particle* part2 = inter.neighbour;
        const Geom::Vector& vec =  inter.distance;
 
        const double pressure2 = part2->get_pressure();
        const double mass2 = part2->get_mass();
        const double rho2  = part2->get_density();
 
        const double h = part2->get_h();
        const Geom::Vector n = Kernel::gk_pres(h, vec);
        force() -=  ((pressure1 + pressure2) / 2.) * (mass2/rho2) * n;
      }
    }
 
    template<class Kernel> void
    Particle::apply_viscous_force() {
      const Geom::Vector& u1 = linear_velocity();
      for (auto& inter: interaction_) {
        SPH::Particle* part2 = inter.neighbour;
        const Geom::Vector& vec =  inter.distance;
 
        const double h = get_h();
        const double mu = part2->get_viscosity();
        const Geom::Vector& u2 = part2->linear_velocity();
        const double mass = part2->get_mass();
        const double rho  = part2->get_density();
        force() += mu*(u2-u1)*(mass/rho)* Kernel::lk_visc(h, vec);
      }
    }
 
    template<class Archive>
    void
    Particle::serialize(Archive& ar, const unsigned int version) {
      ar & boost::serialization::base_object<Physic::Node>(*this);
      ar & _h;
      ar & _radius;
      ar & _density;
      ar & _pressure;
      ar & _viscosity;
    }
 
  }
}
 
 
namespace GranOO3
{
  GRANOO_CLASS_DECLARE_TPL(SPH::Particle);
}
 
 
#endif