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 Diffusion2D {
 
    struct Molecule {
        double diffusionCoef; 
        double defaultQuantity; 
        double defaultEvaporation; 
        Molecule(double dc, double dq, double de)
        : diffusionCoef(dc), defaultQuantity(dq), defaultEvaporation(de) 
        {}
    };
 
    struct GridCell {
        std::vector<double> prevQuantities; 
        std::vector<double> quantities; 
        std::vector<double> consumptions; 
        inline GridCell() = default;
 
        explicit GridCell(std::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:
        std::vector<std::vector<GridCell>> grid; 
        std::vector<Molecule> molecules; 
        bool toreX = true; 
        bool toreY = true; 
        int height = 10;
        int width = 10;
 
        void preStep() {
            for (int i = 0; i < grid.size(); ++i) {
                for (int j = 0; j < grid.size(); ++j) {
                    for(int m = 0; m < molecules.size(); ++m){
                        grid[i][j].prevQuantities[m] = grid[i][j].quantities[m];
                        grid[i][j].quantities[m] = 0.;
                        grid[i][j].consumptions[m] = 0.;
                    }
                }
            }
        }
 
        void diffuse(int x, int y){
            GridCell &gc = grid[x][y];
            int nbNeighbors = 8;
            for (int i = -1; i <= 1; ++i) {
                for (int j = -1; j <= 1; ++j) {
                    if(i != 0 || j != 0){
                        int xi = x + i;
                        int yj = y + j;
                        if (toreX && toreY) { 
                            if(xi < 0) xi = width -1;
                            if(xi >= width) xi = 0;
                            if(yj < 0) yj = height -1;
                            if(yj >= height) yj = 0;
                            for(int m = 0; m < molecules.size(); ++m){
                                gc.quantities[m] += 1./8. * molecules[m].diffusionCoef * grid[xi][yj].prevQuantities[m];
                            }
                        }else if(toreX){
                            if(xi < 0) xi = width -1;
                            if(xi >= width) xi = 0;
                            if(yj >= 0 && yj < height){
                                for(int m = 0; m < molecules.size(); ++m){
                                    gc.quantities[m] += 1./8. * molecules[m].diffusionCoef * grid[xi][yj].prevQuantities[m];
                                }
                            }else nbNeighbors --;
                        }else if(toreY){
                            if(yj < 0) yj = height -1;
                            if(yj >= height) yj = 0;
                            if(xi >= 0 && xi < width){
                                for(int m = 0; m < molecules.size(); ++m){
                                    gc.quantities[m] += 1./8. * molecules[m].diffusionCoef * grid[xi][yj].prevQuantities[m];
                                }
                            }else nbNeighbors --;
                        }else{ 
                            if(xi >= 0 && xi < width && yj >= 0 && yj < height){
                                for(int m = 0; m < molecules.size(); ++m){
                                    gc.quantities[m] += 1./8. * molecules[m].diffusionCoef * grid[xi][yj].prevQuantities[m];
                                }
                            }
                            else nbNeighbors --;
                        }
                    }
                }
            }
            for(int m = 0; m < molecules.size(); ++m){
                gc.quantities[m] += gc.prevQuantities[m] * (1. - molecules[m].diffusionCoef) + (8-nbNeighbors)/8. *molecules[m].diffusionCoef*gc.prevQuantities[m] - molecules[m].defaultEvaporation - gc.consumptions[m];
                if (gc.quantities[m] < 0.) { //can't have a negative amount of molecule
                    gc.quantities[m] = 0.;
                }
            }
        }
 
        void computeStep() {
            for (int i = 0; i < width; ++i) {
                for (int j = 0; j < height; ++j) { //for each cell in the grid
                    diffuse(i,j);
                }
            }
        }
 
        template<typename world_t>
        void updateConsumptions(world_t *w){
            for(auto cell : w->cells){
                int x = cell->getBody().get2DPosition().x();
                int y = cell->getBody().get2DPosition().y();
                for(int m = 0; m < molecules.size(); ++m){
                    grid[x][y].consumptions[m] += cell->getBody().getConsumptions()[m];
                }
            }
        }
 
    public:
 
        std::unordered_map<int,int> moleculesDict;
 
        DiffusionGrid() = default;
 
        inline void setToreX(bool b) { toreX = b; }
 
        inline void setToreY(bool b) { toreY = b; }
 
        inline void setTore(bool x, bool y) { toreX = x; toreY = y; }
 
        inline double getMolecule(int x, int y, int m) { return(grid[x][y].quantities[moleculesDict[m]]); }
 
        inline int getWidth() const { return width; }
 
        inline int getHeight() const { return height; }
 
        inline std::vector<std::vector<GridCell>> &getGrid() { return grid; }
 
        inline std::vector<Molecule> getMolecules() const { return molecules; }
 
        void addMolecule(int id, Molecule m) {
            moleculesDict[id] = (int) molecules.size(); 
            molecules.push_back(m); 
        }
 
        void initGrid(int w, int h){
            grid.resize(w);
            for(int i = 0; i < w; ++i){
                grid[i].resize(h);
                for(int j = 0; j < h; ++j){
                    grid[i][j] = GridCell(molecules);
                }
            }
            height = h;
            width = w;
        }
 
        void initGrid(int size) { initGrid(size, size); }
 
        template<typename world_t>
        void computeMolecules(world_t *w) {
            preStep();
            updateConsumptions(w);
            computeStep();
        }
    };
}
 
#endif //ONKO3D_3_0_DIFFUSIONGRID_HPP