PhysicsConstraints.hpp

Go to the documentation of this file.

#ifndef PHYSICSCONSTRAINTS_HPP
#define PHYSICSCONSTRAINTS_HPP
 
#include <vector>
#include <math.h>
#include <mecacell/mecacell.h>
#include <functional>
 
class PhysicsConstraints {
private:
    class Plane {
    public:
        MecaCell::Vec n; 
        MecaCell::Vec p; 
        bool sticky; 
        Plane(MecaCell::Vec _n = MecaCell::Vec(0, 1, 0), MecaCell::Vec _p = MecaCell::Vec(0, 0, 0), bool _sticky = false) {
            n = _n;
            n.normalize();
            p = _p;
            sticky = _sticky;
        }
 
        double distanceToPoint(MecaCell::Vec c) {
            if (n.x() * c.x() + n.y() * c.y() + n.z() * c.z() + (-n.x() * p.x() - n.y() * p.y() - n.z() * p.z()) < 0)
                return -(MecaCell::Vec::getProjectionOnPlane(p, n, c) - c).length();
            else
                return (MecaCell::Vec::getProjectionOnPlane(p, n, c) - c).length();
        }
    };
    std::vector<Plane> planes;
 
    class SymAxis {
    public:
        MecaCell::Vec axis; 
        MecaCell::Vec origin; 
        double length; 
        double maxEffectRadius; 
        bool fitShape; 
        double attraction_velocity; 
        std::function<double(double, MecaCell::Vec, MecaCell::Vec)> radiusFunction; 
        SymAxis(MecaCell::Vec a, MecaCell::Vec o, double l, double mef, bool fs, double av, std::function<double(double, MecaCell::Vec, MecaCell::Vec)> rf)
            : axis(a), origin(o), length(l), maxEffectRadius(mef), fitShape(fs), attraction_velocity(av), radiusFunction(rf) {}
    };
    std::vector<SymAxis> symAxes;
 
    std::vector<MecaCell::Vec> velocities;
    double ratioPoisson2OnYoung = 0.0000000000001;
    double velocity_center_attraction = 0.00001;
    double radiusBox = -1;
 
public:
    double overlapAdhesion = 0.2;
 
    PhysicsConstraints() {}
 
    void addPlane(MecaCell::Vec _n = MecaCell::Vec(0, 1, 0), MecaCell::Vec _p = MecaCell::Vec(0, 0, 0), bool _sticky = false) {
        planes.push_back(Plane(_n, _p, _sticky));
    }
 
    void addSymAxis(std::function<double(double, MecaCell::Vec, MecaCell::Vec)> _radiusFunction, MecaCell::Vec _axis = MecaCell::Vec(0., 1., 0.), MecaCell::Vec _origin = MecaCell::Vec(0., 0., 0.), double _length = 100, double _maxEffectRadius = 100, bool _fitShape = true, double _velocity_attraction = 0.001) {
        symAxes.push_back(SymAxis(_axis, _origin, _length, _maxEffectRadius, _fitShape, _velocity_attraction, _radiusFunction));
    }
 
    void addVelocity(MecaCell::Vec _v = MecaCell::Vec(0, 1, 0)) {
        velocities.push_back(_v);
    }
 
    template<typename cell_t>
    MecaCell::Vec getSymAxesConstraintsVelocity(cell_t* c) {
        MecaCell::Vec ret = MecaCell::Vec::zero();
        for (auto symAxis : symAxes) {
            double distOnAxis = symAxis.axis.normalized().dot((c->getBody().getPosition() - symAxis.origin));
            // proj dist to origin
            if (distOnAxis > 0 && distOnAxis < symAxis.length) {
                // dist to axis
                MecaCell::Vec dir = (c->getBody().getPosition() - (symAxis.axis.normalized() * distOnAxis + symAxis.origin));
                double sqDistToAxis = dir.sqlength();
                if (sqDistToAxis > 0) dir.normalize();
 
                if (sqDistToAxis < symAxis.maxEffectRadius * symAxis.maxEffectRadius) {
                    double targetRadius = symAxis.radiusFunction(distOnAxis / symAxis.length, symAxis.origin, symAxis.axis.normalized() * symAxis.length);
                    double sqDistToRadius = targetRadius * targetRadius - sqDistToAxis;
                    if ((sqDistToRadius > 0 && symAxis.fitShape) || sqDistToRadius < 0) {
                        // push cell to radius
                        ret = dir * (sqDistToRadius > 0 ? 1 : -1) * symAxis.attraction_velocity;
                    }
                }
            }
        }
        return ret;
    }
 
    template<typename cell_t>
    MecaCell::Vec getConstraintVelocity(cell_t* c) {
        // force pushing cells to center
        double distToCenter = c->getBody().getPosition().length();
        MecaCell::Vec ret = (radiusBox == -1 || distToCenter <= radiusBox)
            ? MecaCell::Vec::zero()
            : -c->getBody().getPosition().normalized() * std::pow(distToCenter - radiusBox, 2.0) * velocity_center_attraction;
        // constraints velocities
        for (auto v : velocities)
            ret += v;
        // shape constraints
        ret += getSymAxesConstraintsVelocity(c);
        return ret;
    }
 
    void setCurvedBox(double radiusBoxPlane, double curvedDistToMaxSpeed, double maxSpeed) {
        radiusBox = radiusBoxPlane;
        velocity_center_attraction = maxSpeed / std::pow(curvedDistToMaxSpeed, 2.0);
    }
 
    void setPoissonAndYoungOfPlane(double poisson, double young) {
        ratioPoisson2OnYoung = (1.0 - poisson * poisson) / young;
    }
 
    template<typename cell_t>
    void applyPlaneConstraints(cell_t* c) {
        double dist, delta, R, E, Frep, adhesion, Fadh;
        for (auto plane : planes) {
            if (plane.sticky) {
                dist = std::abs(plane.distanceToPoint(c->getBody().getPosition()));
                if (dist <= c->getBody().getBoundingBoxRadius()) {
                    delta = c->getBody().getBoundingBoxRadius() - dist;
                    R = 1.0 / (1.0 / c->getBody().getBoundingBoxRadius() + 1.0 / c->getBody().getBoundingBoxRadius());
                    E = 1.0 /
                        ((1.0 - c->getBody().getPoissonNumber() * c->getBody().getPoissonNumber()) / c->getBody().getYoungModulus() +
                        ratioPoisson2OnYoung);
 
                    Frep = 4.0 / 3.0 * E * sqrt(R) * pow(delta, 3.0 / 2.0);
                    adhesion = (4.0 * sqrt(overlapAdhesion * c->getBody().getBoundingBoxRadius()) * E * sqrt(R)) / (3.0 * M_PI * M_PI * R * R);
                    Fadh = -M_PI * M_PI * R * R * adhesion * delta;
 
                    c->getBody().receiveForce(plane.n * (Frep + Fadh));
                }
            }
        }
    }
 
    template<typename cell_t>
    void snapCellsToPlane(cell_t* c) {
        double dist;
        for (auto plane : planes) {
            dist = plane.distanceToPoint(c->getBody().getPosition()) - (c->getBody().getBoundingBoxRadius() - c->getBody().getBoundingBoxRadius() * overlapAdhesion);
            if (dist < 0) c->getBody().setPosition(c->getBody().getPosition() - plane.n * dist);
        }
    }
};
 
#endif // PHYSICSCONSTRAINTS_HPP