
// c0 lap u + c1 u = f 
// CL u|x=0 = 0
// f = g
// time scheme backward euler 
// [ c0 masse + c1 rigidity ] u^n = fv + rho*epsilon/tau masse  V^n-1
// + [ rho*epsilon/dt^2*(1 + dt alpha) masse + beta*\lambda1/tau rigidity ] u^{n-1}

 
// parameters
real density = 0.9, young = 0.9e5, poisson = 0.4, thickness = 0.1;
real alpha = 0., beta = 0; //(coeff rayleigh)
real dt = 0.0001, iMax = 500.;
//real tgv = 1e30; // a hude value of exact penalisation of boundary condition

func h= sin(pi*x);

real lambda0 = young*thickness/(0.25*(1-poisson*poisson)); 
real lambda1 = young*thickness/(2.*(1.+poisson));   
real c0 = density*thickness/(dt*dt)+lambda0+alpha*density*thickness/dt;
real c1 = lambda1+beta*lambda1/dt;
real a = density*thickness*(1.+dt*alpha)/(dt*dt);
real b = density*thickness/dt;

macro Grad(u) [dx(u),dy(u),dz(u)]  // EOM

// string mesh 
meshS Th=square3(20,50,[x,3.*y,z]);
savemesh(Th,"ThL.mesh");
// FE space
fespace Uh(Th,P1);
Uh u, v, V, uold ;
V=0;uold=0;

func g=10000.; //x^2;//; //10000.
// c0 * mass matrix + c1 stifness matrix
varf m(u,v) = int2d(Th) (u*v);
varf k(u,v) = int2d(Th) (Grad(u)'*Grad(v));
varf mmkk(u,v) = int2d(Th) (c0*u*v + c1*Grad(u)'*Grad(v)) + on(2,u=0);
varf cl(u,v) = on(2,u=0);

matrix mass=m(Uh,Uh);
matrix stifness=k(Uh,Uh);
matrix mk = mmkk(Uh,Uh); //,tgv=tgv,solver="SPARSESOLVER");

varf f(u,v) = int2d(Th)(g*v);
real[int] Force=f(0,Uh);
real[int] CL=cl(0,Uh);
matrix T = a*mass + (beta*lambda1/dt)*stifness;  

//initial condition
u[] = 0;
  
meshS Thmv;
real t=dt;

for (int i=0; i<=iMax; i++) {
	   V[]=(u[]-uold[]);
	   V[]/=dt; 
	   uold[]=u[];
	   //  for the time dependant part
	   real[int] rhs = mass*V[];  	   
	   rhs *= b;
	   rhs += CL;
	   rhs += T*uold[];
	   // apply the external force
	   if(i==1 || i==100)
	   		rhs+=Force;
	   u[] = mk^-1*rhs;
	   t+=dt;
	   // coeff for the visu
	   real coeff=100;
	   Thmv=movemesh(Th,[x,y,u*coeff]);
	   if (verbosity)
	   		cout << " time " << t << " u min "<< u[].min << " u max " << u[].max << endl;
       plot(Th, Thmv, wait=0, cmm=" d = "+t+" iter = "+i,prev=1);


    }

