//dvx/dt =dsxx/dx+dsxy/dy (momentum conservation in x)

//dvy/dt =dsxy/dx+dsyy/dy (momentum conservation in y)

//Structure:

//dsxx/dt=dvx/dx; dsxy/dt=(dvx/dy+dvy/dx)/2, dsyy/dt=dvy/dy

//Fluid:(nu viscosity + nonlinear terms in momentum eq)

////sxx= nu*dvx/dx; sxy=nu*(dvx/dy+dvy/dx)/2, //syy=nu*dvy

let vx = [];

let vy = [];

let mx = [];

let my = [];

let x = [];

let y = []; 

let sxx = [];

let syy= [];

let sxy= [];

let N=100;

let t=0;

let dt;

let h;

let nu;

let S=1;

​

​

function setup() {

  createCanvas(400, 400);

  for (var i=0;i<N;i++){

    x[i]=[];

    y[i]=[];

    vx[i]=[];

    vy[i]=[];

    sxx[i]=[];

    syy[i]=[];

    sxy[i]=[];

  for (var j=0;j<N;j++){

    x[i][j]=4*i;

    y[i][j]=4*j;

    vx[i][j]=0

    vy[i][j]=0;

    sxx[i][j]=0;

    syy[i][j]=0;

    sxy[i][j]=0;

  }

  }

  for (var i=35;i<65;i++){

  for (var j=45;j<55;j++){

  vy[i][j]=10;

  vx[i][j]=0;

  }

  }

  h=1/N;

  nu=0.2*h

  dt=0.1*h;

}

​

function draw() {

  background(255,255,255,100);

//  text("Vector Burgers: 2d Shock Problem",100,100);

  noStroke();

  t=t+dt;

//  fill(0);

  for (var i=1;i<N-1;i++){

  for (var j=1;j<N-1;j++){

//  if (j<N/2){

//  sxx[i][j]=sxx[i][j]+dt*(vx[i][j]-vx[i-1][j])/h;

//  sxy[i][j]=sxy[i][j]+dt*0.5*(vx[i][j]-vx[i][j-1]+vy[i][j]-vy[i-1][j])/h;

// syy[i][j]=syy[i][j]+dt*(vy[i][j]-vy[i][j-1])/h;

//  vx[i][j]=vx[i][j]+0.1*dt*(vx[i+1][j]+vx[i-1][j]-4*vx[i][j]+vx[i][j+1]+vx[i][j-1])/h+dt*(sxx[i+1][j]-sxx[i][j]+sxy[i][j+1]-sxy[i][j])/h;

//  vy[i][j]=vy[i][j]+0.1*dt*(vy[i+1][j]+vy[i-1][j]-4*vy[i][j]+vy[i][j+1]+vy[i][j-1])/h+dt*(sxy[i+1][j]-sxy[i][j]+syy[i][j+1]-syy[i][j])/h;

// }

//  else{

  sxx[i][j]=nu*(vx[i][j]-vx[i-1][j])/h;

  sxy[i][j]=nu*0.5*(vx[i][j]-vx[i][j-1]+vy[i][j]-vy[i-1][j])/h;

  syy[i][j]=nu*(vy[i][j]-vy[i][j-1])/h;

  vx[i][j]=vx[i][j]+0.5*dt*vx[i][j]*(vx[i+1][j]-vx[i-1][j])/h+0.5*dt*vy[i][j]*(vx[i][j+1]-vx[i][j-1])/h+1*dt*(vx[i+1][j]+vx[i-1][j]-4*vx[i][j]+vx[i][j+1]+vx[i][j-1])/h+dt*(sxx[i+1][j]-sxx[i][j]+sxy[i][j+1]-sxy[i][j])/h;

vy[i][j]=vy[i][j]+0.5*dt*vx[i][j]*(vy[i+1][j]-vy[i-1][j])/h+0.5*dt*vy[i][j]*(vy[i][j+1]-vy[i][j-1])/h+1*dt*(vy[i+1][j]+vy[i-1][j]-4*vy[i][j]+vy[i][j+1]+vy[i][j-1])/h+dt*(sxy[i+1][j]-sxy[i][j]+syy[i][j+1]-syy[i][j])/h;

}

}

​

fill(0);

if (t>0.1){

for (var i=0;i<N;i++){

for (var j=0;j<N;j++){

  x[i][j]=4*i;

  y[i][j]=4*j;

}

}

t=0;

}

​

for (var i=0;i<N;i++){

for (var j=0;j<N;j++){

  x[i][j]=x[i][j]+200*dt*vx[i][j];

  y[i][j]=y[i][j]+200*dt*vy[i][j];

  if (j<N/2){

  fill(255,0,0,50*(pow(vx[i][j],2)+pow(vy[i][j],2)));

  ellipse(4*i,4*j,3);

  }

  if (j>N/2){

  fill(0,0,255,255);

  ellipse(x[i][j],y[i][j],2);

  }

  }

  }

 fill(0);

for (var i=1;i<N-1;i++){

fill(0);

ellipse(4*i,200-20*(pow(vx[i][45],2)+pow(vy[i][45],2)),3);

fill(0,255,0,255);

ellipse(4*i,y[i][50],3);

}

fill(0);

}