DiffusionGrid.hpp

Go to the documentation of this file.

#ifndef ONKO3D_3_0_DIFFUSIONGRID_HPP
#define ONKO3D_3_0_DIFFUSIONGRID_HPP
 
#include <vector>
#include <typeindex>
#include <typeinfo>
#include <unordered_map>
#include <mecacell/utilities/utils.hpp>
 
namespace Diffusion {
 
    struct Molecule {
        static constexpr double spheroidDensity = 0.001; 
        static constexpr double henryConst = 22779000000.; 
        double diffusion; 
        double defaultQuantity; 
        double omega; 
        double defaultEvaporation; 
        inline Molecule(double d, double dQ, double density, double dE) : diffusion(d), defaultQuantity(dQ),
        omega(spheroidDensity  * density * henryConst), defaultEvaporation(dE) {}
    };
 
    struct GridCell {
        vector<double> prevQuantities; 
        vector<double> quantities; 
        vector<double> consumptions; 
        inline GridCell() = default;
 
        explicit GridCell(vector<Molecule> m) {
            int size = m.size();
            quantities.resize(size);
            prevQuantities.resize(size);
            consumptions.assign(size, 0.);
            for (int i = 0; i < size; ++i) {
                quantities[i] = m[i].defaultQuantity;
                prevQuantities[i] = m[i].defaultQuantity;
            }
        }
    };
 
    class DiffusionGrid {
    private:
        double dx = 5.; 
        double accuracy = 0.1; 
        unordered_map<MecaCell::Vec, GridCell> grid; 
        vector<Molecule> molecules; 
        // Boundaries
        MecaCell::Vec minBoundary;
        MecaCell::Vec maxBoundary;
 
        inline double deriv2nd(double xm1, double x, double xp1) { return ((xm1 - (2.0 * x) + xp1) / (dx * dx)); }
 
        void nextStep(int i) {
            for (auto &c : grid) {
                c.second.prevQuantities[i] = c.second.quantities[i];
            }
        }
 
        void allCellEmpty() {
            for (auto &c : grid) {
                for (int i = 0; i < molecules.size(); ++i) { // For each molecule
                    c.second.consumptions[i] = 0.;
                }
            }
        }
 
        void cleanCellOutBoundary() {
            for (auto &c: grid) {
                if (isOutBoundary(c.first)) {
                    grid.erase(c.first);
                }
            }
        }
 
        bool isOutBoundary(const MecaCell::Vec &v) {
            return(v.x()<minBoundary.x() || v.y()<minBoundary.y() || v.z()<minBoundary.z() || v.x()>maxBoundary.x() || v.y()>maxBoundary.y() || v.z()>maxBoundary.z());
        }
 
        inline double getDt(double D) const { return (dx * dx) / (6 * D); }
 
        double computeStep(int i) {
            Molecule m = molecules[i]; // Considered molecule
            double D = m.diffusion; // Diffusion constant in µm²/s
            double dt = getDt(D); // Diffusion time step in s
            double total = 0; // Total amount of considered molecule in mmHg
            for (auto &c : grid) { // For each cell in the grid
                GridCell &gc = c.second;
                double prevQuantities = gc.prevQuantities[i];
                // Laplacian in mmHg/µm²
                double laplacian = deriv2nd(getMolecule(c.first - MecaCell::Vec(1, 0, 0), i), prevQuantities,
                                            getMolecule(c.first + MecaCell::Vec(1, 0, 0), i)) +
                                   deriv2nd(getMolecule(c.first - MecaCell::Vec(0, 1, 0), i), prevQuantities,
                                            getMolecule(c.first + MecaCell::Vec(0, 1, 0), i)) +
                                   deriv2nd(getMolecule(c.first - MecaCell::Vec(0, 0, 1), i), prevQuantities,
                                            getMolecule(c.first + MecaCell::Vec(0, 0, 1), i));
                gc.quantities[i] = prevQuantities + (D * laplacian - ((gc.consumptions[i] + m.defaultEvaporation) * m.omega)) * dt;
                if (gc.quantities[i] < 0.) { // Can't have a negative amount of molecule
                    gc.quantities[i] = 0.;
                }
                total += gc.quantities[i];
            }
            return (total);
        }
 
    public:
 
        unordered_map<int,int> moleculesDict; 
        DiffusionGrid() = default;
 
        inline DiffusionGrid(double dx, double accuracy) : dx(dx), accuracy(accuracy), grid(), molecules() {}
 
        inline double getDx() const { return dx; }
 
        inline void setDx(double _dx) { dx = _dx; }
 
        inline void setAccuracy(double a) { accuracy = a; }
 
        MecaCell::Vec getIndexFromPosition(const MecaCell::Vec &v) const {
            MecaCell::Vec res = v / dx;
            return MecaCell::Vec(floor(res.x()), floor(res.y()), floor(res.z()));
        }
 
        double getMolecule(MecaCell::Vec v, int m) {
            if (grid.count(v) <= 0)
                return (molecules[m].defaultQuantity);
            else
                return(grid[forward<MecaCell::Vec>(v)].quantities[m]);
        }
 
        inline size_t size() const { return grid.size(); }
 
        inline unordered_map<MecaCell::Vec, GridCell> &getGrid() { return grid; }
 
        inline double getCellSize() const { return dx; }
 
        inline vector<Molecule> getMolecules() const { return molecules; }
 
        void addMolecule(int n, Molecule m) {
            moleculesDict[n] = (int) molecules.size(); 
            molecules.push_back(m); 
        }
 
        template<typename world_t>
        void updateBoundary(world_t *w) {
            MecaCell::Vec min(9999,9999,9999);
            MecaCell::Vec max(-9999,-9999,-9999);
            allCellEmpty();
            for (auto &c : w->cells) { // For each cell
                MecaCell::Vec center = getIndexFromPosition(c->getBody().getPosition());
 
                // Update min and max boundaries
                if(center.x()<min.x()) min.setX(center.x());
                if(center.y()<min.y()) min.setY(center.y());
                if(center.z()<min.z()) min.setZ(center.z());
 
                if(center.x()>max.x()) max.setX(center.x());
                if(center.y()>max.y()) max.setY(center.y());
                if(center.z()>max.z()) max.setZ(center.z());
                
                // Update consumptions and create new grid cells
                if (grid.count(forward<MecaCell::Vec>(center)) <= 0) {
                    GridCell gc = GridCell(molecules);
                    for (int i = 0; i < molecules.size(); ++i) { // For each molecule
                        gc.consumptions[i] += c->getBody().getConsumptionByIndex(i);
                    }
                    grid[forward<MecaCell::Vec>(center)] = gc;
                } else {
                    for (int i = 0; i < molecules.size(); ++i) { // For each molecule
                        grid[forward<MecaCell::Vec>(center)].consumptions[i] += c->getBody().getConsumptionByIndex(i);
                    }
                }
            }
 
            minBoundary = min;
            maxBoundary = max;
            for(int x=minBoundary.x();x<maxBoundary.x();x++){
                for(int y=minBoundary.y();y<maxBoundary.y();y++){    
                    for(int z=minBoundary.z();z<maxBoundary.z();z++){
                        MecaCell::Vec pos(x,y,z);
                        if (grid.count(forward<MecaCell::Vec>(pos)) <= 0) {
                            grid[forward<MecaCell::Vec>(pos)] = GridCell(molecules);
                        }
                    }   
                }
            }
        }
 
        template<typename world_t>
        void computeMolecules(world_t *w) {
            vector<double> lasts; // Last step quantities
            vector<double> news; // Current step quantities
            lasts.assign(molecules.size(),0.);
            news.assign(molecules.size(), INFINITY);
            for (int i = 0; i < molecules.size(); ++i) { // For each molecule
                double t = 0; // Elapsed time
                double dt = getDt(molecules[i].diffusion); // Diffusion time step
                while (abs(lasts[i] - news[i]) / grid.size() >= accuracy && t < w->dt) {
                    lasts[i] = news[i];
                    nextStep(i);
                    news[i] = computeStep(i);
                    t += dt;
                }
            }
        }
 
        inline double getMoleculeRealPos(const MecaCell::Vec &v, int mol) {
            return (getMolecule(getIndexFromPosition(v), mol));
        }
    };
}
 
#endif //ONKO3D_3_0_DIFFUSIONGRID_HPP