BodyHertzianPhysics.hpp

Go to the documentation of this file.

#ifndef BODYMASSSPRINGDAMPER_HPP
#define BODYMASSSPRINGDAMPER_HPP
#define MOVABLE
 
#include <mecacell/mecacell.h>
#include <mecacell/movable.h>
#include "PluginHertzianPhysics.hpp"
#include "../../../../src/core/BaseBody.hpp"
#include <math.h>
#include "Tensor.hpp"
 
namespace HertzianPhysics {
 
    template<typename cell_t, class plugin_t>
    class BodyHertzianPhysics : public MecaCell::Movable, public virtual BaseBody<plugin_t> {
 
    private:
        double baseRadius; 
        double radius; 
        double volume; 
        double density; 
        double gammaCoef; 
        MecaCell::Vec purposeVelocity; 
        double poissonNumber; 
        double youngModulus; 
        double adhesion; 
        double targetOverlap; 
        bool computed;
        Tensor stressTensor; 
        inline void computeAdhesion() {
            double E = 1.0 / (2.0 * (1.0 - poissonNumber * poissonNumber) / youngModulus);
            double R = 1.0 / (2.0 / radius);
            adhesion = (4.0 * sqrt(targetOverlap * radius) * E * sqrt(R)) /
                    (3.0 * M_PI * M_PI * R * R);
        }
 
    public:
        
        BodyHertzianPhysics()
        : stressTensor()
        {
            baseRadius = 10.0; // 10 µm radius
            radius = baseRadius;
            volume = 4.0 * M_PI / 3.0 * pow(radius, 3); // µm³
            density = 0.001; // Density of water ~ 1000 kg/m³ = 10¹⁵ ng/m³ = 10^-3 ng/µm³ 
            poissonNumber = 0.5; // Coefficient between 0 and 1
            youngModulus = 4000000000.0; // 4 kPa = 10³ kg/(m.s²) = 10¹⁵ ng/(m.s²) = 10⁹ ng/(µm.s²) 
            this->setMass(density * volume); // Mass = density * volume µg
            gammaCoef = 100000000000.0;
            purposeVelocity = MecaCell::Vec::zero();
            setTargetOverlap(0.2);
        }
 
        inline double getBoundingBoxRadius() const { return radius; }
 
        inline void setRadius(double rad) {
            radius = rad;
            volume = 4.0 * M_PI / 3.0 * pow(radius, 3.0);
            this->setMass(density * volume); // Mass = density * 4π/3 * radius³ ng
            computeAdhesion();
        }
 
        inline double getVolume() const { return volume; }
 
        inline double getRadius() const { return radius; }
 
        inline double getBaseRadius() const { return baseRadius; }
 
        inline void setBaseRadius(double _baseRadius) { baseRadius = _baseRadius; }
 
        inline void growth(double delta) {
            this->setRadius(baseRadius * std::cbrt(1 + delta));
        }
 
        inline void setDensity(double d) {
            density = d;
            this->setMass(density * volume); // Mass = density * 4π/3 * radius³ ng
        }
 
        inline double getAdhesion() const { return adhesion; }
 
        inline void setTargetOverlap(double _targetOverlap) {
            this->targetOverlap = _targetOverlap;
            computeAdhesion();
        }
 
        inline double getYoungModulus() const { return youngModulus; }
 
        inline void setYoungModulus(double ym) { youngModulus = ym; }
 
        inline MecaCell::Vec getPurposeVelocity() const { return purposeVelocity; }
 
        inline void setPurposeVelocity(MecaCell::Vec pf) { purposeVelocity = pf; }
 
        inline double getPoissonNumber() const { return poissonNumber; }
 
        inline void setPoissonNumber(double pn) { poissonNumber = pn; }
 
        inline void moveTo(const MecaCell::Vec &v) { this->setPosition(v); }
 
        inline bool getComputed() const { return computed; }
 
        inline void setComputed(bool c) { computed = c; }
 
        inline double getGamma() const { return this->getMass() * gammaCoef; }
 
        inline void setGammaCoef(double g) { gammaCoef = g; }
 
        inline Tensor getStressTensor() { return stressTensor; }
 
        inline void setStressTensor(Tensor st) { stressTensor = st; }
 
        inline void addStress(Tensor st) { stressTensor += st; }
 
        inline void divTensor(double div) { stressTensor /= div; }
 
        inline double getPressure() const { return abs((stressTensor.get(0, 0) + stressTensor.get(1, 1) + stressTensor.get(2, 2)) / 3.0); }
    };
}
 
#endif // BODYMASSSPRINGDAMPER_HPP