// To generate movie set movie=true
// You should have ImageMagick and ffmepg installed

bool movie=true;

verbosity=0;
real L=8,H=3;			// half length and height of the domain
real r=1,gamma=0.5;		// diameter of the ellipse, anisotropy
real theta=pi/8.,s1=-5,s2=0;	// theta=angle of rotation of the solid. s=[s1,s2]=translation of the solid
real mu=5.e-3,rho=1;		// viscosity and density of the fluid
real rho0=10;			// density of the solid
macro g [0,0]			// gravitational potential
real dt=0.1,T=100;		// time step, Final time

real uD=(x-H)^2/H^2;
cout<<"Reynolds="<<rho*uD*L/mu<<endl;

// definition of the mesh //
int Dirichlet=1,Gamma=2,inflow=3,outflow=4;
border bottom(t=-L,L){x=t;y=-H;label=Dirichlet;}
border right(t=-H,H){x=L;y=t;label=outflow;}
border top(t=L,-L){x=t;y=H;label=Dirichlet;}
border left(t=H,-H){x=-L;y=t;label=inflow;}
border ball(t=0,2.*pi){x=r*(cos(theta)*cos(t)-gamma*sin(theta)*sin(t))+s1;y=r*(gamma*cos(theta)*sin(t)+sin(theta)*cos(t))+s2;label=Gamma;}
real dn=10,h=L/dn;	// scale of the mesh
mesh Th=buildmesh(bottom(dn*2*L)+right(dn*2*H)+top(dn*2*L)+left(2*dn*H)+ball(2*r*pi*dn));
int Ball=Th(s1,s2).region;	// label of the region that contains the ball
int Fluid=Th(L,H).region;	// label of the fluid region

real I=int2d(Th,Ball)(x'*x);	// Inertia
real V=int2d(Th,Ball)(1.);	// Volume


fespace Vh(Th,[P2,P2]),Qh(Th,P1),Rh(Th,P1dc);	// FE spaces for the velocity/pressure/rotational
real delta=1.e-5,epsilon=0;		// Penalization of e(u)=0 on the Ball and of the pression

// Some useful macros
macro s [s1,s2]//
macro u [u1,u2]//
macro up [up1,up2]//
macro v [v1,v2]//
real sqrt2=sqrt(2.);
macro e(u) [dx(u[0]),(dx(u[1])+dy(u[0]))/sqrt2,dy(u[1])]//
macro div(u) (dx(u[0])+dy(u[1]))//
macro rot(u) (dx(u[1])-dy(u[0]))//
macro wedge(u,v) ((u)[0]*(v)[1]-(u)[1]*(v)[0])//
macro Convect(u,dt,v) [convect((u),(dt),(v)[0]),convect((u),(dt),(v)[1])]//

Vh u,up,v; Qh p,q;
real t;
problem Stokes(u1,u2,p,v1,v2,q,init=t)=
	 int2d(Th,Fluid)(rho/dt*u'*v)
	-int2d(Th,Fluid)(rho/dt*Convect(up,-dt,up)'*v)
	
	+int2d(Th,Ball)(rho0/dt*u'*v)
	-int2d(Th,Ball)(rho0/dt*up'*v)
	
	+int2d(Th)(mu*e(u)'*e(v) + div(u)*q+div(v)*p - epsilon*p*q)
	+int2d(Th,Ball)(e(u)'*e(v)/delta)
	-int2d(Th,Fluid)(rho*g'*v)
	-int2d(Th,Ball)(rho0*g'*v)

+on(Dirichlet,u1=0,u2=0)
+on(inflow,u1=uD,u2=0)
+on(Gamma,u1=0,u2=0);

real[int] viso(50); for(int i=0;i<viso.n;i++) viso(i)=-70+140/(viso.n-1)*i;
real[int] colors=[
1, 1 , 0,  //
1, 1 , 0.5,  //
1, 1 , 1.,  //
0.5, 1 , 1,  //
0.5, 1 , 0.5,  //
0.5, 1 , 0.,  //
0.5, 1 , 0.5,  //
0.5, 1 , 1,  //
0., 1 , 1.,  //
0., 1 , 0.5,  //
0., 1 , 0];  //

up=[uD,0];	// initialization
string input="convert ";
int i=0;
if (movie)
for(t=0;t<T;t+=dt)
{
	//Th=buildmesh(bottom(dn*2*L)+right(dn*2*H)+top(dn*2*L)+left(2*dn*H)+ball(2*r*pi*dn));
	Ball=Th(s1,s2).region;
	Stokes;
	Rh omega=rot(u);
	if(!movie) plot(omega,cmm="time t="+t+"min="+omega[].min+";max="+omega[].max,fill=1,viso=viso,hsv=colors);
	if (movie)
	{
		plot(omega,fill=1,viso=viso,hsv=colors,ps="tmp/rotational-"+i+".eps");
		//file<<"tmp/rotational-"+i+".eps ";
		input=input+"tmp/rotational-"+i+".eps ";
	}
	i++;
	up=u;
}
//file<<"tmp/NS.jpg"<<endl;
string inter="tmp/NS.jpg";
if (movie)
{
	exec("convert "+input+inter);
	exec("ffmpeg -i tmp/NS-%d.jpg NS.mpg");
	exec("ffmpeg -i tmp/NS-%d.jpg NS.ogg");
	exec("ffmpeg -i tmp/NS-%d.jpg NS.webm");
	exec("ffmpeg -i tmp/NS-%d.jpg NS.swf");
	exec("ffmpeg -i tmp/NS-%d.jpg -vcodec libx264 -vf 'scale=854:-2' -b 854k  NS.mp4");
//	exec("rm -rf tmp/*");
}