> show canvas only <


/* built with Studio Sketchpad: 
 *   https://sketchpad.cc
 * 
 * observe the evolution of this sketch: 
 *   https://studio.sketchpad.cc/sp/pad/view/ro.B0W-rseERXA/rev.4553
 * 
 * authors: 
 *   
 *   haha sait nimen
 *   
 *   atomim
 *   

 * license (unless otherwise specified): 
 *   creative commons attribution-share alike 3.0 license.
 *   https://creativecommons.org/licenses/by-sa/3.0/ 
 */ 



// Pressing Control-R will render this sketch.

int i = 0; 

Simulation simulation;

void setup() {  // this is run once.   
    simulation = new Simulation();
    // set the background color
    background(255);
    
    // canvas size (Variable aren't evaluated. Integers only, please.)
    size(300, 300);

    // limit the number of frames per second
    frameRate(30);
} 

void draw() {  // this is run repeatedly.  
    simulation.update();
    simulation.draw();
}

class GridPoint{
    PVector velocity;
    PVector newVelocity;
    float mass;
    float denisity;
    boolean hasParticles;
    GridPoint(){
        this.velocity=new PVector(0,0);
        this.newVelocity=new PVector(0,0);
    }
}

class Grid{
    PVector gravity;
    GridPoint[][] grid;
    int w;
    int h;
    Grid(int w, int h){
        this.gravity=new PVector(0,0.0003);
        this.grid=new GridPoint[w][h];
        for(x=0;x<w;x++){
            for(y=0;y<h;y++){
                this.grid[x][y]=new GridPoint();
            }
        }
        this.w=w;
        this.h=h;
        
    }
    
    void draw(){
        
        for (int y=0;y<this.h;y++){
            for (int x=0;x<this.w;x++){
                var stepx=(width/w);
                var stepy=(height/h);
                var point=this.grid[x][y];
                fill(point.mass*200,point.velocity.x*2000+128, point.velocity.y*2000+128);
                rect(x*stepx,y*stepy,stepx,stepy);
                
            }
        }
    }
    void reset(){
        for (int y=0;y<this.h;y++){
            for (int x=0;x<this.w;x++){
                var point=this.grid[x][y];
                point.mass=0;
                point.velocity.x=0;point.velocity.y=0;
            }
        }
    }
    
    float totalMass(){
        float sum=0;
        for (int y=0;y<this.h;y++){
            for (int x=0;x<this.w;x++){
                var point=this.grid[x][y];
                sum+=point.mass;
            }
        }
        return sum;
    }
    
    void updateVelocities(){
        //just gravity
        for (int y=0;y<this.h;y++){
            for (int x=0;x<this.w;x++){
                this.grid[x][y].newVelocity.y=this.grid[x][y].velocity.y+this.gravity.y;
                this.grid[x][y].newVelocity.x=this.grid[x][y].velocity.x;
            }
        }
    }
}

class Particle{
    PVector velocity;
    PVector position;
    Particle(x,y){
        this.velocity=new PVector(random(-0.04,0.04),random(-0.04,0.04));
        this.position=new PVector(x,y);
    }
}

class Particles{
    PVector gravity;
    Particle[] particles;
    int amount;
    Grid grid;
    Particles(Grid grid){
        this.particles=new Particle[0];
        this.amount=0;
        this.gravity=new PVector(0,0.0003);
        this.grid=grid;
    }
    void draw(){
        fill(255);
        for(int i=0;i<this.particles.length;i++){
            Particle particle=this.particles[i];
        rect(particle.position.x*width/grid.w, particle.position.y*height/grid.h,5,5)
        }
    }
    void add(particle){
        println(this.particles.length);
        this.particles[this.particles.length++]=particle
    }
    
    
    void interpolateMassAndVelocityTo(Grid grid){
        //suppose weight support size is one
        //suppose particle mass equals 1
        for(int i=0;i<this.particles.length;i++){
            Particle particle=this.particles[i];
            for(int x=max(0,floor(particle.position.x-0.5)); x<=min(ceil(particle.position.x-0.5),grid.w-1); x++){
                for(int y=max(0,floor(particle.position.y-0.5)); y<=min(ceil(particle.position.y-0.5),grid.h-1); y++){
                    grid.grid[x][y].mass+=weight(particle.position,x,y);
                }
            }
        }
        for(int i=0;i<this.particles.length;i++){
            Particle particle=this.particles[i];
            for(int x=max(0,floor(particle.position.x-0.5)); x<=min(ceil(particle.position.x-0.5),grid.w-1); x++){
                for(int y=max(0,floor(particle.position.y-0.5)); y<=min(ceil(particle.position.y-0.5),grid.h-1); y++){
                    int wht=weight(particle.position,x,y)/grid.grid[x][y].mass; //normalized to conserve momentum
                    grid.grid[x][y].velocity.x+=particle.velocity.x*wht;
                    grid.grid[x][y].velocity.y+=particle.velocity.y*wht;
                }
            }
        }
    }
    
    void updateVelocitiesFrom(grid){
        float picAmount=0.05;
        //PIC/FLIP update
        for(int i=0;i<this.particles.length;i++){
            Particle particle=this.particles[i];
            float flipx=0;
            float flipy=0;
            float picx=0;
            float picy=0;
            for(int x=max(0,floor(particle.position.x-0.5)); x<=min(ceil(particle.position.x-0.5),grid.w-1); x++){
                for(int y=max(0,floor(particle.position.y-0.5)); y<=min(ceil(particle.position.y-0.5),grid.h-1); y++){
                    int wht=weight(particle.position,x,y);
                    flipx+=(grid.grid[x][y].newVelocity.x-grid.grid[x][y].velocity.x)*wht
                    flipy+=(grid.grid[x][y].newVelocity.y-grid.grid[x][y].velocity.y)*wht
                    picx+=(grid.grid[x][y].newVelocity.x)*wht
                    picy+=(grid.grid[x][y].newVelocity.y)*wht
                }
            }
            particle.velocity.x=(particle.velocity.x+flipx)*(1-picAmount)+picAmount*picx;
            particle.velocity.y=(particle.velocity.y+flipy)*(1-picAmount)+picAmount*picy;
        }
    }
    
    void advance(){
        
        for(int i=0;i<this.particles.length;i++){
            this.particles[i].position.add(this.particles[i].velocity);
        }
    }
}

class Simulation {
    Particles particles;
    Grid grid;
    
    Simulation(){
        this.grid= new Grid(10,10);
        this.particles=new Particles(grid);
    }
    
    void draw(){
        grid.draw();
        particles.draw();
    }
    
    
    // https://www.math.ucla.edu/~jteran/papers/SSCTS13.pdf
    void update(){
        grid.reset();
        if(int(random(10))==1)this.particles.add(new Particle(grid.w/2+random(-0.5,0.5),1.5+random(-0.5,0.5)));

        // 1) complete using simplified interpolation
        particles.interpolateMassAndVelocityTo(grid);
        // 2) Not relevant yet
        // 3)
        // 4) just gravity
        grid.updateVelocities();
        // 8) complete PIC/FLIP
        particles.updateVelocitiesFrom(grid);
        // 10)
        particles.advance();

        println(grid.totalMass());
    }
}


float weight(PVector position, int x, int y){
    //linear interpolation to start with.
    //gridposition is supposed to be integers
    return max(0,1-0.5*abs((x-position.x+0.5))-0.5*abs((y-position.y+0.5)))*0.5;
}

/*
TODO: 
-add timestep size
*/